[R] survreg penalized likelihood?

[Thomas Lumley]

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