R Under development (unstable) (2025-10-28 r88973 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.

> 
> library("mlt")
Loading required package: basefun
Loading required package: variables
> library("survival")
> set.seed(29)
> options(digits = 3)
> 
> 
> ### true dgp
> rY <- function(n, ...) rexp(n, ...)
> pY <- function(x, ...) pexp(x, ...)
> dY <- function(x, ...) dexp(x, ...)
> 
> ### three groups
> gf <- gl(3, 1)
> g <- rep(gf, 500)
> 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.556  1.220  1.220  1.220  1.220  1.898  0.839  1.222 
> coef(cph <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata))
   g2    g3 
0.693 1.004 
> 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.0392  0.2923  0.6962  1.0579  1.0579  1.7642  0.6928  1.0042 
> ## 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.7480 -0.4207  0.0848  0.4320  0.4320  1.0904  0.6935  1.0066 
> coef(cph <- coxph(y ~ g, data = mydata))
   g2    g3 
0.693 1.007 
> 
> ### 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.413  0.630  1.177  1.177  1.333  1.945  0.540  0.732 
> 
> ### 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.2736 -0.0105  0.7839  0.7839  0.9045  1.4596  0.5717  0.6999 
> 
> ### uncensored, time-varying coefficients in both groups
> mydata <- data.frame(y = y, g = g)
> m <- ctm(response = logBb, 
+            interacting = as.basis(~ 0 + g, data = mydata),
+            todist = "MinExtrVal")
> logLik(opt <- mlt(m, data = mydata, scale = TRUE))
'log Lik.' -682 (df=18)
> logLik(cph <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata))
'log Lik.' -9359 (df=2)
> ## 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 
   2.46    0.32    2.78