[R] survreg penalized likelihood?
Spencer Graves
spencer.graves at pdf.com
Mon Apr 21 01:33:44 CEST 2003
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
>
>
More information about the R-help
mailing list