R Under development (unstable) (2023-06-29 r84618 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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. > > library("mlt") Loading required package: basefun Loading required package: variables > library("survival") > set.seed(29) > options(digits = 5) > > > ### true dgp > rY <- function(n, ...) rexp(n, ...) > pY <- function(x, ...) pexp(x, ...) > dY <- function(x, ...) dexp(x, ...) > > ### tree groups > gf <- gl(3, 1) > g <- rep(gf, 100) > y <- rY(length(g), rate = (1:nlevels(g))[g]) > mydata <- data.frame(y = y, g = g) > > boxplot(y ~ g, data = mydata) > > ### uncensored, Cox model, h = bernstein > Bb <- Bernstein_basis(numeric_var("y", support = c(0, max(y) + .1), bounds = c(0, Inf)), + order = 5, ui = "increasing") > s <- as.basis(~ g, data = data.frame(g = gf), remove_intercept = TRUE) > m <- ctm(response = Bb, shifting = s, todist = "MinExtrVal") > (cf1 <- coef(opt <- mlt(m, data = mydata))) Bs1(y) Bs2(y) Bs3(y) Bs4(y) Bs5(y) Bs6(y) g2 g3 -2.57065 1.10359 1.10359 1.10359 1.10359 1.80043 0.93238 1.26275 > coef(cph <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata)) g2 g3 0.69981 1.02270 > yn <- mkgrid(Bb, 50)$y > yn <- yn[yn > 0] > a <- predict(opt, newdata = data.frame(g = gf[1]), q = yn) > layout(matrix(1:4, ncol = 2)) > plot(yn, a, type = "l", col = "red") > lines(yn, log(yn)) > a <- predict(opt, newdata = data.frame(g = gf), q = yn, type = "survivor") > plot(yn, a[,1], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[1]))) > plot(yn, a[,2], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[2]))) > plot(yn, a[,3], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[3]))) > > ### h = bernstein(log()) > logBb <- Bernstein_basis(numeric_var("y", support = c(1, max(y) + .1), bounds = c(min(y) / 2, Inf)), + order = 5, ui = "increasing", log_first = TRUE) > m <- ctm(response = logBb, shifting = s, todist = "MinExtrVal") > (cf1 <- coef(opt <- mlt(m, data = mydata))) Bs1(y) Bs2(y) Bs3(y) Bs4(y) Bs5(y) Bs6(y) g2 g3 -0.11781 0.19444 0.75901 0.90914 0.90914 1.65722 0.74255 1.03796 > ## sample from this model > sam <- simulate(opt, newdata = data.frame(g = gf), nsim = 100) > nd <- data.frame(y = unlist(sam), g = rep(gf, length(sam))) > opt2 <- mlt(m, data = nd) > ## visualise > yn <- mkgrid(Bb, 50)$y > yn <- yn[yn > 0] > a <- predict(opt, newdata = data.frame(g = gf[1]), q = yn) > layout(matrix(1:4, ncol = 2)) > plot(yn, a, type = "l", col = "red") > lines(yn, log(yn)) > a <- predict(opt, newdata = data.frame(g = gf), q = yn, type = "survivor") > plot(yn, a[,1], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[1]))) > plot(yn, a[,2], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[2]))) > plot(yn, a[,3], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[3]))) > > ### right censoring > mydata <- data.frame(y = Surv(y, sample(0:1, length(y), replace = TRUE)), g = g) > coef(opt <- mlt(m, data = mydata, scale = TRUE)) Bs1(y) Bs2(y) Bs3(y) Bs4(y) Bs5(y) Bs6(y) g2 g3 -0.80017 -0.50675 -0.30211 0.25437 0.25437 1.20222 0.66696 1.19417 > coef(cph <- coxph(y ~ g, data = mydata)) g2 g3 0.59631 1.15907 > > ### left censoring > mydata <- data.frame(y = Surv(y, sample(0:1, length(y), replace = TRUE), type = "left"), g = g) > coef(opt <- mlt(m, data = mydata, scale = TRUE)) Bs1(y) Bs2(y) Bs3(y) Bs4(y) Bs5(y) Bs6(y) g2 g3 0.43810 0.64728 1.05050 1.05050 1.05050 1.79448 0.48891 0.66483 > > ### interval censoring > mydata <- data.frame(y = Surv(y, y + 1, sample(0:3, length(y), replace = TRUE), type = "interval"), + g = g) > coef(opt <- mlt(m, data = mydata, scale = TRUE)) Bs1(y) Bs2(y) Bs3(y) Bs4(y) Bs5(y) Bs6(y) g2 g3 -0.48006 -0.18524 0.26866 0.57464 0.57464 1.27249 0.80125 0.98864 > > ### uncensored, time-varying coefficients in both groups > mydata <- data.frame(y = y, g = g) > m <- ctm(response = logBb, + interacting = as.basis(~ g, data = mydata), + todist = "MinExtrVal") > ## IGNORE_RDIFF_BEGIN > op <- mltoptim(spg = list(maxit = 5000, quiet = TRUE, checkGrad = FALSE)) > coef(opt <- mlt(m, data = mydata, optim = op, scale = TRUE)) Bs1(y):(Intercept) Bs2(y):(Intercept) Bs3(y):(Intercept) Bs4(y):(Intercept) -0.12013 0.14785 0.78542 0.84945 Bs5(y):(Intercept) Bs6(y):(Intercept) Bs1(y):g2 Bs2(y):g2 0.84945 1.63768 0.74501 0.79693 Bs3(y):g2 Bs4(y):g2 Bs5(y):g2 Bs6(y):g2 0.78117 0.71714 0.71714 2.67694 Bs1(y):g3 Bs2(y):g3 Bs3(y):g3 Bs4(y):g3 1.08890 1.15693 0.51936 1.70777 Bs5(y):g3 Bs6(y):g3 1.70777 30.41933 > ## IGNORE_RDIFF_END > coef(cph <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata)) g2 g3 0.69981 1.02270 > ## visualize > a <- predict(opt, newdata = data.frame(g = gf[1]), q = yn) > layout(matrix(1:4, ncol = 2)) > plot(yn, a, type = "l", col = "red") > lines(yn, log(yn)) > a <- predict(opt, newdata = data.frame(g = gf), q = yn, type = "survivor") > plot(yn, a[,1], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[1]))) > plot(yn, a[,2], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[2]))) > plot(yn, a[,3], type = "l", col = "red", ylim = c(0, 1)) > lines(survfit(cph, newdata = data.frame(g = gf[3]))) > > > proc.time() user system elapsed 3.07 0.59 3.65