Estimating shape and scale parameters of Weibull distribution, with right censoring

[Franklin Dexter]

I have a continuous variable to which I want to fit a Weibull distribution. The variable is right censored. My objective is to estimate the shape and scale parameters. Because I’m failing with my data, I created a hypothetical example:

set obs 5000
generate x = rweibull(6.183, 12.23)

generate xCensor = 0
replace xCensor = 1 if runiform() < 0.01

stset x, failure(xCensor)
streg x, distribution(weibull)



_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x | .0069157 .0043201 -7.96 0.000 .0020329 .0235266
_cons | 9.62e-39 1.00e-37 -8.41 0.000 1.34e-47 6.93e-30
-------------+----------------------------------------------------------------
/ln_p | 4.054029 .1238296 32.74 0.000 3.811327 4.29673
-------------+----------------------------------------------------------------
p | 57.62916 7.136195 45.2104 73.4592
1/p | .0173523 .0021487 .013613 .0221188
------------------------------------------------------------------------------

The estimated shape parameter 4.054 differs considerably from 6.183. In addition, the estimated scale parameter, exp(0.0069157) = 1.

I’d appreciate assistance in understanding what I’m clearly doing incorrectly to estimate the 6.183 and 12.23.

[Paul Dickman]

Nobody else has responded, so I'll give a couple of brief comments. I don't have Stata on this machine so can't provide code.

You are generating uncensored survival times with the given distribution. You are then modifying the data, so that for 90% of the observations the true survival time is no longer the survival time but a right censored observation of the survival time. That is, the true value of the survival time will be greater than the observed (censored) value. As such, after you've changed the data they are no longer a sample from the original data generating process.

When fitting the Weibull model, you should not put x in the model as a predictor. x is the outcome (the survival time).

Not an error, but the code is easier to read if xCensor takes the value 1 for censored observations; you have specified that observations with xCensor==1 are failures.

There are lots of good papers showing how to simulate survival data. e.g. https://www.stata-journal.com/articl...article=st0275 or https://journals.sagepub.com/doi/pdf...867X1201200405

[Franklin Dexter]

My objective is not the simulation, I have data. I'm using the simulation only to show the problem. I revised accordingly and increased sample size:

set obs 50000
generate x = rweibull(6.183, 12.23)

generate xCensor = 1
replace xCensor = 0 if runiform() < 0.01

stset x, failure(xCensor)
streg, distribution(weibull)


------------------------------------------------------------------------------
_t | Hazard Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 1.86e-07 1.04e-08 -277.53 0.000 1.67e-07 2.08e-07
-------------+----------------------------------------------------------------
/ln_p | 1.821697 .0034979 520.79 0.000 1.814842 1.828553
-------------+----------------------------------------------------------------
p | 6.182344 .0216255 6.140104 6.224875
1/p | .161751 .0005658 .1606458 .1628637
------------------------------------------------------------------------------

My objective is not the simulation, just to estimate shape 6.183 and scale 12.23. This is giving the 6.183. How obtain scale 12.23?

[Paul Lambert]

There are different ways to define parameters in distributions and it is always important to check how any software has implemented things.
If you use rweibullph() rather than rweibull() you will get the same parameterization as used in streg.

[Paul Dickman]

I think the following code illustrates what you want.

set obs 50000
set seed 123456
generate x = rweibullph(6.183, 12.23)

generate xCensor = 1

stset x, failure(xCensor) exit(time 0.7)
streg, distribution(weibull)

I added a seed so you'll get the exact same output.

Output:

------------------------------------------------------------------------------
          _t |     Hazard   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   12.31991   .1573629   196.60   0.000     12.01531    12.63223
-------------+----------------------------------------------------------------
       /ln_p |   1.823543   .0045531   400.51   0.000     1.814619    1.832467
-------------+----------------------------------------------------------------
           p |   6.193765   .0282005                      6.138739    6.249285
         1/p |   .1614527   .0007351                      .1600183    .1628999
------------------------------------------------------------------------------

Using the exit() option to -stset- is a better way of enforcing censoring. The mean of the distribution is around 0.6 so I set all survival times greater than 0.7 to be censored. That is, uncensored survival times of, for example, 0.8 are replaced with 0.7 and the failure indicator set to 0. This resulted in 25.8 percent of the survival times being censored.

. tab _d

       1 if |
 failure; 0 |
if censored |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |     12,914       25.83       25.83
          1 |     37,086       74.17      100.00
------------+-----------------------------------
      Total |     50,000      100.00