{smcl}
{com}{sf}{ul off}{txt}{.-}
       log:  {res}C:\Documents and Settings\jli2\My Documents\Log\MLE_log\duration.smcl
  {txt}log type:  {res}smcl
 {txt}opened on:  {res} 2 Oct 2004, 16:55:31
{txt}
{com}. * duration analysis_Weibull model
. * use the data of stata official website
. clear
{txt}
{com}. use http://www.stata-press.com/data/r8/kva
{txt}(Generator experiment)

{com}. * canned routine
. * the default model is the Weibull model with proportional hazard
. list
{txt}
     {c TLC}{hline 10}{c -}{hline 6}{c -}{hline 10}{c TRC}
     {c |} {res}failtime   load   bearings {txt}{c |}
     {c LT}{hline 10}{c -}{hline 6}{c -}{hline 10}{c RT}
  1. {c |} {res}     100     15          0 {txt}{c |}
  2. {c |} {res}     140     15          1 {txt}{c |}
  3. {c |} {res}      97     20          0 {txt}{c |}
  4. {c |} {res}     122     20          1 {txt}{c |}
  5. {c |} {res}      84     25          0 {txt}{c |}
     {c LT}{hline 10}{c -}{hline 6}{c -}{hline 10}{c RT}
  6. {c |} {res}     100     25          1 {txt}{c |}
  7. {c |} {res}      54     30          0 {txt}{c |}
  8. {c |} {res}      52     30          1 {txt}{c |}
  9. {c |} {res}      40     35          0 {txt}{c |}
 10. {c |} {res}      55     35          1 {txt}{c |}
     {c LT}{hline 10}{c -}{hline 6}{c -}{hline 10}{c RT}
 11. {c |} {res}      22     40          0 {txt}{c |}
 12. {c |} {res}      30     40          1 {txt}{c |}
     {c BLC}{hline 10}{c -}{hline 6}{c -}{hline 10}{c BRC}

{com}. * we hope to fit a Weibull model in which failtime depends on load and bearings
. * first claim failtime as [st] data
. stset failtime

     {txt}failure event:  {res}(assumed to fail at time=failtime)
{txt}obs. time interval:  {res}(0, failtime]
{txt} exit on or before:  {res}failure

{txt}{hline 78}
{res}       12{txt}  total obs.
{res}        0{txt}  exclusions
{hline 78}
{res}       12{txt}  obs. remaining, representing
{res}       12{txt}  failures in single record/single failure data
{res}      896{txt}  total analysis time at risk, at risk from t = {res}        0
                             {txt}earliest observed entry t = {res}        0
                                  {txt}last observed exit t = {res}      140
{txt}
{com}. * run the regression
. streg load bearings, distribution(weibull)

         {txt}failure _d:  {res}1{txt} (meaning all fail)
   analysis time _t:  {res}failtime

{txt}Fitting constant-only model:

Iteration 0:   log likelihood = {res}-13.666193
{txt}Iteration 1:   log likelihood = {res}-9.7427276
{txt}Iteration 2:   log likelihood = {res}-9.4421169
{txt}Iteration 3:   log likelihood = {res}-9.4408287
{txt}Iteration 4:   log likelihood = {res}-9.4408286

{txt}Fitting full model:

Iteration 0:   log likelihood = {res}-9.4408286{txt}  
Iteration 1:   log likelihood = {res} -2.078323{txt}  
Iteration 2:   log likelihood = {res} 5.2226016{txt}  
Iteration 3:   log likelihood = {res} 5.6745808{txt}  
Iteration 4:   log likelihood = {res} 5.6934031{txt}  
Iteration 5:   log likelihood = {res} 5.6934189{txt}  
Iteration 6:   log likelihood = {res} 5.6934189{txt}  

Weibull regression -- log relative-hazard form 

No. of subjects = {res}          12                     {txt}Number of obs   ={res}        12
{txt}No. of failures = {res}          12
{txt}Time at risk    = {res}         896
                                                   {txt}LR chi2({res}2{txt})      ={res}     30.27
{txt}Log likelihood  =   {res} 5.6934189                     {txt}Prob > chi2     ={res}    0.0000

{txt}{hline 13}{c TT}{hline 64}
          _t {c |} Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
{hline 13}{c +}{hline 64}
        load {c |}  {res} 1.599315   .1883807     3.99   0.000     1.269616    2.014631
    {txt}bearings {c |}  {res} .1887995   .1312109    -2.40   0.016     .0483546    .7371644
{txt}{hline 13}{c +}{hline 64}
       /ln_p {c |}  {res} 2.051552   .2317074     8.85   0.000{col 58} 1.597414    2.505691
{txt}{hline 13}{c +}{hline 64}
           p {c |}  {res} 7.779969   1.802677{col 58} 4.940241    12.25202
         {txt}1/p {c |}  {res} .1285352   .0297826{col 58} .0816192    .2024193
{txt}{hline 13}{c BT}{hline 64}

{com}. * the results above only reports hazard ratio, but we need coefficients. Hence we add nohr option
. streg load bearings, distribution(weibull) nohr

         {txt}failure _d:  {res}1{txt} (meaning all fail)
   analysis time _t:  {res}failtime

{txt}Fitting constant-only model:

Iteration 0:   log likelihood = {res}-13.666193
{txt}Iteration 1:   log likelihood = {res}-9.7427276
{txt}Iteration 2:   log likelihood = {res}-9.4421169
{txt}Iteration 3:   log likelihood = {res}-9.4408287
{txt}Iteration 4:   log likelihood = {res}-9.4408286

{txt}Fitting full model:

Iteration 0:   log likelihood = {res}-9.4408286{txt}  
Iteration 1:   log likelihood = {res} -2.078323{txt}  
Iteration 2:   log likelihood = {res} 5.2226016{txt}  
Iteration 3:   log likelihood = {res} 5.6745808{txt}  
Iteration 4:   log likelihood = {res} 5.6934031{txt}  
Iteration 5:   log likelihood = {res} 5.6934189{txt}  
Iteration 6:   log likelihood = {res} 5.6934189{txt}  

Weibull regression -- log relative-hazard form 

No. of subjects = {res}          12                     {txt}Number of obs   ={res}        12
{txt}No. of failures = {res}          12
{txt}Time at risk    = {res}         896
                                                   {txt}LR chi2({res}2{txt})      ={res}     30.27
{txt}Log likelihood  =   {res} 5.6934189                     {txt}Prob > chi2     ={res}    0.0000

{txt}{hline 13}{c TT}{hline 64}
          _t {c |}      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
{hline 13}{c +}{hline 64}
        load {c |}  {res} .4695753   .1177884     3.99   0.000     .2387143    .7004363
    {txt}bearings {c |}  {res}-1.667069   .6949745    -2.40   0.016    -3.029194   -.3049443
       {txt}_cons {c |}  {res}-45.13191   10.60663    -4.26   0.000    -65.92053   -24.34329
{txt}{hline 13}{c +}{hline 64}
       /ln_p {c |}  {res} 2.051552   .2317074     8.85   0.000{col 58} 1.597414    2.505691
{txt}{hline 13}{c +}{hline 64}
           p {c |}  {res} 7.779969   1.802677{col 58} 4.940241    12.25202
         {txt}1/p {c |}  {res} .1285352   .0297826{col 58} .0816192    .2024193
{txt}{hline 13}{c BT}{hline 64}

{com}. * user-defined MLE procedure
. program drop _all
{txt}
{com}. program define duration_stata
{txt}  1{com}. args lnf theta1 theta2
{txt}  2{com}. quietly replace `lnf'=`theta1'+ln(`theta2')+(`theta2'-1)*ln($ML_y1)-exp(`theta1')*$ML_y1^`theta2'
{txt}  3{com}. end
{txt}
{com}. ml model lf duration_stata ( failtime= load bearings) ()
{txt}
{com}. ml maximize

{txt}initial:       log likelihood = {res}    -<inf>{txt}  (could not be evaluated)
feasible:      log likelihood = {res}-192.82771
{txt}rescale:       log likelihood = {res}-90.249353
{txt}rescale eq:    log likelihood = {res}-78.855988
{txt}Iteration 0:   log likelihood = {res}-78.855988{txt}  
Iteration 1:   log likelihood = {res}-56.436521{txt}  
Iteration 2:   log likelihood = {res}-48.035487{txt}  
Iteration 3:   log likelihood = {res}-47.523484{txt}  (not concave)
Iteration 4:   log likelihood = {res}-44.634008{txt}  
Iteration 5:   log likelihood = {res}-44.401497{txt}  
Iteration 6:   log likelihood = {res}-44.396795{txt}  
Iteration 7:   log likelihood = {res}-44.396793{txt}  

{col 51}Number of obs{col 67}= {res}        12
{col 51}{txt}Wald chi2({res}2{txt}){col 67}= {res}     15.98
{txt}Log likelihood = {res}-44.396793{col 51}{txt}Prob > chi2{col 67}= {res}    0.0003

{txt}{hline 13}{c TT}{hline 64}
    failtime {c |}      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
{hline 13}{c +}{hline 64}
{res}eq1          {txt}{c |}
        load {c |}  {res} .4695755   .1177995     3.99   0.000     .2386928    .7004583
    {txt}bearings {c |}  {res} -1.66707   .6949983    -2.40   0.016    -3.029242   -.3048986
       {txt}_cons {c |}  {res}-45.13193   10.60782    -4.25   0.000    -65.92287     -24.341
{txt}{hline 13}{c +}{hline 64}
{res}eq2          {txt}{c |}
       _cons {c |}  {res} 7.779973   1.802881     4.32   0.000     4.246392    11.31355
{txt}{hline 13}{c BT}{hline 64}

{com}. * note we have the same results as canned routine.
. log close
       {txt}log:  {res}C:\Documents and Settings\jli2\My Documents\Log\MLE_log\duration.smcl
  {txt}log type:  {res}smcl
 {txt}closed on:  {res} 2 Oct 2004, 17:00:17
{txt}{.-}
{smcl}
{txt}{sf}{ul off}