[R] survreg penalized likelihood?

[Spencer Graves]

Thanks very much.

After I distributed that example, I noticed that the documtation for the 
survival package says, "Survival analysis, including penalised 
likelihood", so I thought that some penalty function might be invoked 
under certain circumstances.

Thanks again,
Spencer Graves

Thomas Lumley wrote:
> On Sat, 19 Apr 2003, Spencer Graves wrote:
> 
> 
>>	  What objective function is maximized by survreg with the default
>>Weibull model?  I'm getting finite parameters in a case that has the
>>likelihood maximzed at Infinite, so it can't be a simple maximum
>>likelihood.
> 
> 
> The objective function *is* the loglikelihood.
> 
> 
>>Consider the following:
>>#############################
>> > set.seed(3)
>> > Stress <- rep(1:3, each=3)
>> > ch.life <- exp(9-3*Stress)
>> > simLife <- rexp(9, rate=1/ch.life)
>> > Data <- data.frame(Stress=factor(Stress),
>>+  Time=pmin(simLife, 50), dead = (simLife<=50))
>> > Data[1:3,]
>>   Stress Time  dead
>>1      1   50 FALSE
>>2      1   50 FALSE
>>3      1   50 FALSE
>># All three observations at Stress=1 are censored.
>> > Fit <- survreg(Surv(Time, dead)~Stress-1, data=Data)
>> > Fit
>>Call:
>>survreg(formula = Surv(Time, dead) ~ Stress - 1, data = Data)
>>
>>Coefficients:
>>    Stress1    Stress2    Stress3
>>41.3724178  2.2819204 -0.5281005
>>
>>Scale= 0.2577886
>>
>>Loglik(model)= -6.6   Loglik(intercept only)= -22.1
>>	Chisq= 30.88 on 2 degrees of freedom, p= 2e-07
>>n= 9
>># Even though 100% of observations at Stress=1 are censored,
>># I still get a finite estimate for log(characteristics life).
>>
> 
> 
> 
> I suppose technically 41.3 is a finite estimate of log(life), but since
> your data are censored at 50, a life expectancy of nearly 10^18 isn't
> terribly finite.  The likelihood at this value is very very very close to
> the likelihood at the true mle, and that's all a numerical optimisation
> technique can really be expected to give you (and all that statistical
> theory, frequentist or Bayesian, demands).
> 
> It's quite difficult to come up with a reliable test for ridges in the
> loglikelihood -- coxph() tries, but gives too many false positives.
> 
> Incidentally, if you relax the assumption that the scale parameter is the
> same for each Stress you end up with much larger finite predictions for
> Stress=1, as the scale parameter ends up around 10^15 in that stratum.
> 
> 	-thomas
> 
>