[R] Fitting of lognormal distribution to lower tail experimental data
Göran Broström
goran.brostrom at gmail.com
Sat Jan 17 00:31:17 CET 2009
On Sat, Jan 17, 2009 at 12:27 AM, Göran Broström
<goran.brostrom at gmail.com> wrote:
> On Fri, Jan 16, 2009 at 3:31 PM, Mattias Brännström
> <Mattias.Brannstrom at tt.luth.se> wrote:
>> Thank you, David!
>>
>> I agree and apprechiate your analysis, which definitely will influence my
>> analysis of this data, but still I would like you to disregard from it(!)
>>
>> The standard routine in the field is, beyond my control, to assume
>> lognormal distribution to achieve comparable results also with other
>> materials (comparison is made on COV). For that reason I have to use it,
>> even if it is not statistically defendable for this particular data.
>>
>> So, if I rephrase the question to be (more general):
>> How would you fit a lognormal distribution to the lower 10% tail of the
>> data (assuming it was lognormal)? What functions to use?
>
> Mattias, it is not clear (to me) what you mean by "fit a lognormal
> distribution to the 10%-lower tail of the attached data" (and what is
> COV?). However, a guess is that you really mean what you say, so I
> tried to right-censor your data at the 10% quantile (33.4134, Type I
> censoring)
Sorry, it is of course censoring of Type II.
Göran
and fit the resulting data to a lognormal distribution. The
> fit was fairly good, as can be seen by comparing the fitted cumulative
> hazard function to the corresponding non-parametric one (the
> Nelson-Aalen estimator):
>
>> cc <- read.table(....
>> v1 <- ifelse(cc$V1 <= 33.4134, cc$V1, 33.4134)
>> event <- as.numeric(cc$V1 <= 33.4134)
>> library(eha)
> Loading required package: survival
> Loading required package: splines
>> fit.ln <- phreg(Surv(v1, event) ~ 1, dist = "lognormal")
>> fit.cox <- coxreg(Surv(v1, event) ~ 1)
>> check.dist(fit.cox, fit.ln)# Gives you a plot
>> summary(fit.ln)
> Call:
> phreg(formula = Surv(v1, event) ~ 1, dist = "lognormal")
>
> Covariate Coef Exp(Coef) se(Coef) Wald p
> log(scale) 3.943 51.597 0.053 0.000
> log(shape) 1.089 2.970 0.101 0.000
>
> Events 70
> Total time at risk 22998
> Max. log. likelihood -401.38
>
> Is this maybe what you are looking for?
>
> HTH (!)
>
> Göran
>
>> Best regards,
>> Mattias
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
>
>
> --
> Göran Broström
>
--
Göran Broström
More information about the R-help
mailing list