R Under development (unstable) (2025-01-06 r87534 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> require("fitdistrplus")
Loading required package: fitdistrplus
Loading required package: MASS
Loading required package: survival
>
> compareplotdistcens <- function(d)
+ {
+ par(mfrow = c(2,2))
+ plotdistcens(d, NPMLE.method = "Turnbull.middlepoints")
+ plotdistcens(d, NPMLE.method = "Turnbull.intervals")
+ plotdistcens(d, NPMLE.method = "Wang")
+
+ }
>
> comparenpmle <- function(d)
+ {
+ npmleT <- fitdistrplus:::npmle(d, method = "Turnbull.intervals")
+ npmleW <- fitdistrplus:::npmle(d, method = "Wang")
+ print(npmleT)
+ print(npmleW)
+ cat("nb of intervals for Turnbull: ", nrow(npmleT), "\n")
+ cat("nb of intervals for Wang: ", nrow(npmleW), "\n")
+ xmin <- min(c(npmleT$left[is.finite(npmleT$left)],npmleT$left[is.finite(npmleT$left)]))
+ xmax <- max(c(npmleT$right[is.finite(npmleT$right)],npmleT$right[is.finite(npmleT$right)]))
+ par(mfrow = c(2, 1))
+ plot(npmleT$left, npmleT$p, xlim = c(xmin, xmax), main = "left")
+ points(npmleW$left, npmleW$p, col = "red", pch = 4)
+ plot(npmleT$right, npmleT$p, xlim = c(xmin, xmax), main = "right")
+ points(npmleW$right, npmleW$p, col = "red", pch = 4)
+ }
> #### Comparison of plotdistcens with different NPMLE methods
> if(FALSE)
+ {
+
+ # d1 = trivial case with only interval censored data
+ d1 <- data.frame(left = c(1, 2, 3, 4, 3, 7), right = c(2, 5, 3, 7, 8, 9))
+ d <- d1
+ par(mfrow = c(2,2))
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "norm")
+ fb <- fitdistcens(d, "logis")
+ par(mfrow = c(2,2))
+ cdfcompcens(list(fa, fb), NPMLE.method = "Turnbull.middlepoints")
+ cdfcompcens(list(fa, fb), NPMLE.method = "Turnbull.intervals")
+ cdfcompcens(list(fa, fb), NPMLE.method = "Wang")
+ par(mfrow = c(1,2))
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+
+ # d2 = case with left and right censored data
+ data(smokedfish)
+ d2 <- smokedfish
+ d <- d2
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "lnorm")
+ fb <- fitdistcens(d, "gamma")
+ par(mfrow = c(2,2))
+ cdfcompcens(list(fa, fb), NPMLE.method = "Turnbull.middlepoints")
+ cdfcompcens(list(fa, fb), NPMLE.method = "Turnbull.intervals")
+ cdfcompcens(list(fa, fb), NPMLE.method = "Wang")
+ par(mfrow = c(1,2))
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+ # d3 = case with also rigth censored data
+ d3 <- data.frame(left = c(-1.4, 1.18, -1.4, 2, -1.4, 0),
+ right = c(1, 1.18, 2, NA, 0, 2))
+ d <- d3
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "norm")
+ fb <- fitdistcens(d, "logis")
+ par(mfrow = c(1,2))
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+ # d4 = case with also right censored data
+ # with differences between the algorithms by the way they polish
+ # the ECDF function, by putting some masses to zero.
+ require("actuar")
+ data(fluazinam)
+ d4 <- fluazinam
+ d <- d4
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "lnorm")
+ fb <- fitdistcens(d, "llogis")
+ par(mfrow = c(1,2))
+ cdfcompcens(list(fa, fb), NPMLE.method = "Turnbull.intervals")
+ cdfcompcens(list(fa, fb), NPMLE.method = "Wang")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+
+ # d5 a random example with exact values
+ set.seed(123)
+ r <- rnorm(10)
+ d5 <- data.frame(left = r, right = r)
+ d <- d5
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "norm")
+ fb <- fitdistcens(d, "logis")
+ par(mfrow = c(1,2))
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+ # d7 = bigger dataset with also rigth censored data
+ data(salinity)
+ d7 <- log10(salinity)
+ d <- d7
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "logis")
+ fb <- fitdistcens(d, "norm")
+ par(mfrow = c(1,2))
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+
+ # d8 = an random example with all types of data (a small one)
+ # set.seed(1234) # check OK
+ # set.seed(1231) # check OK
+ set.seed(1232)
+ ns <- 25
+ r <- rnorm(ns)
+ d8 <- data.frame(left = r, right = r)
+ delta <- rlnorm(ns)
+ icensored <- rbinom(ns, size = 1, prob = 0.2)
+ Lcensored <- rbinom(ns, size = 1, prob = 0.2*(1 - icensored))
+ Rcensored <- rbinom(ns, size = 1, prob = 0.3*(1 - icensored)*(1 - Lcensored))
+ # icensored + Lcensored + Rcensored
+ d8$left <- d8$left * (1 - Lcensored) + (-1000) * Lcensored
+ d8$right <- d8$right * (1 - Rcensored) + (1000) * Rcensored
+ d8$right <- d8$right + delta * icensored
+ d8$right[d8$right == 1000] <- NA
+ d8$left[d8$left == -1000] <- NA
+ d8
+ d <- d8
+ compareplotdistcens(d)
+ comparenpmle(d)
+ fa <- fitdistcens(d, "logis")
+ fb <- fitdistcens(d, "norm")
+ par(mfrow = c(1,2))
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ ppcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Turnbull.intervals")
+ qqcompcens(list(fa, fb), ynoise = FALSE, NPMLE.method = "Wang")
+
+ # d8 = an random example with all types of data (a big one)
+ # set.seed(1234) # check OK
+ # set.seed(1231) # check OK
+ set.seed(1232)
+ ns <- 500
+ # ns <- 5000
+ r <- rnorm(ns)
+ d8 <- data.frame(left = r, right = r)
+ delta <- rlnorm(ns)
+ icensored <- rbinom(ns, size = 1, prob = 0.2)
+ Lcensored <- rbinom(ns, size = 1, prob = 0.2*(1 - icensored))
+ Rcensored <- rbinom(ns, size = 1, prob = 0.3*(1 - icensored)*(1 - Lcensored))
+ # icensored + Lcensored + Rcensored
+ d8$left <- d8$left * (1 - Lcensored) + (-1000) * Lcensored
+ d8$right <- d8$right * (1 - Rcensored) + (1000) * Rcensored
+ d8$right <- d8$right + delta * icensored
+ d8$right[d8$right == 1000] <- NA
+ d8$left[d8$left == -1000] <- NA
+ d <- d8
+
+ par(mfrow = c(2,2))
+ system.time(plotdistcens(d, NPMLE.method = "Turnbull.middlepoints"))
+ system.time(plotdistcens(d, NPMLE.method = "Turnbull.intervals"))
+ system.time(plotdistcens(d, NPMLE.method = "Wang"))
+
+ }
>
> proc.time()
user system elapsed
1.14 0.10 1.23