In R, how do I plot a Weibull probability density function for right censored data using package fitdistrplus?

231 Views Asked by At

I am trying to overlay a Weibull probabiliy density function (PDF) for right censored data on a histogram of the data using package fitdistrplus but have been unable to do so.

df is a data.frame containing a dummy data set in the format of my larger data set.

> library(fitdistrplus)
> left     <- c(rep(1,12),rep(2,5),rep(3,3),rep(4,2),rep(5,1))
> right    <- c(rep(1,12),rep(2,5),rep(3,3),rep(4,2),NA)
> (df      <- data.frame(left,right))
   left right
1     1     1
2     1     1
3     1     1
4     1     1
5     1     1
6     1     1
7     1     1
8     1     1
9     1     1
10    1     1
11    1     1
12    1     1
13    2     2
14    2     2
15    2     2
16    2     2
17    2     2
18    3     3
19    3     3
20    3     3
21    4     4
22    4     4
23    5    NA
> fitcen   <- fitdistrplus::fitdistcens(df,"weibull")
> fitdistrplus::plotdistcens(df,
+                            distr = "weibull",
+                            para = list(fitcen$estimate[1],fitcen$estimate[2]))

Executing the script above generates the figures below. enter image description here

I have successfully overlaid a Weibull PDF over a histogram of a non-censored version of this dummy data set (below). The top-left figure in the image below is very nearly my target except that it shows the PDF of non-censored rather than right censored data.

> noncens <- fitdistrplus::fitdist(left, distr = "weibull")
> plot(noncens)

enter image description here

1

There are 1 best solutions below

0
overcup On

Edit: as I read more I see that this is likely a multi-step process and more complex than a single answer.

Partial answer: extract slope and scale values from fitcen and use use function dweibull() to generate plot of PDF. Histogram of values is not here included.

> shape <- fitcen$estimate[1]
> scale <- fitcen$estimate[2]
> curve(dweibull(x, shape=shape, scale=scale), from=0, to=10)

enter image description here