R Under development (unstable) (2025-10-07 r88904 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("eha")
> library("survival")
> set.seed(29)
> options(digits = 4)
> 
> ## ************** Exponential - AFT *********************************
> 
> ### true dgp
> rY <- function(n, ...) rexp(n, ...)
> pY <- function(x, ...) pexp(x, ...)
> dY <- function(x, ...) dexp(x, ...)
> 
> 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)
> 
> Bb <- log_basis(numeric_var("y", support = range(y)), ui = "increasing",
+                 remove_intercept = TRUE)
> Bx <- as.basis(~ g, data = mydata)
> m <- ctm(Bb, shifting = Bx, todist = "MinExtrVal")
> 
> ## Estimate coefficients
> coef(opt <- mlt(m, data = mydata, fixed = c("log(y)" = 1)))
     log(y) (Intercept)          g2          g3 
     1.0000     -0.1230      0.7915      1.1191 
> 
> coef(aft <- survreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata,
+                     dist = "exponential"))
(Intercept)          g2          g3 
     0.1230     -0.7915     -1.1191 
> 
> coef(cox <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata))
    g2     g3 
0.6998 1.0227 
> 
> coef(phreg <- phreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata, 
+                     dist = "weibull", shape = 1))
        g2         g3 log(scale) 
    0.7915     1.1191     0.1230 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1414 0.1414 
> ## Cox
> sqrt(diag(vcov(cox)))
    g2     g3 
0.1523 0.1557 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1414 0.1414 
> 
> ## Compare log-Likelihoods
> logLik(aft)
'log Lik.' -145.8 (df=3)
> logLik(opt)
'log Lik.' -145.8 (df=3)
> 
> 
> ## Use a Weibull-AFT for estimation (Weibull shape parameter should be nu = 1)
> 
> ## Estimate coefficients
> (cf <- coef(opt2 <- mlt(m, data = mydata)))
     log(y) (Intercept)          g2          g3 
    0.92384    -0.07597     0.72340     1.02745 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
   -0.08223     0.78303     1.11215 
> 
> coef(aft2 <- survreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata,
+                      dist = "weibull"))
(Intercept)          g2          g3 
    0.08223    -0.78304    -1.11215 
> 
> ## Compare Weibull shape paramters
> 1 / cf[1]
log(y) 
 1.082 
> aft2$scale
[1] 1.082
> 
> ## Compare log-Likelihoods
> logLik(opt2)
'log Lik.' -144.3 (df=4)
> logLik(aft2)
'log Lik.' -144.3 (df=4)
> 
> sqrt(diag(vcov(opt2)))[c("g2", "g3")]
    g2     g3 
0.1463 0.1502 
> sqrt(diag(vcov(aft2)))[c("g2", "g3")]
    g2     g3 
0.1532 0.1531 
> 
> 
> ## *************** Right-censored
> 
> mydata <- data.frame(y = Surv(y, sample(0:1, length(y), replace = TRUE)), g = g)
> coef(opt <- mlt(m, data = mydata, fixed = c("log(y)" = 1)))
     log(y) (Intercept)          g2          g3 
     1.0000     -0.7579      0.8466      1.2263 
> 
> ## Estimate coefficients 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
     0.7579     -0.8466     -1.2263 
> coef(cox <- coxph(y ~ g, data = mydata))
    g2     g3 
0.8333 1.2134 
> coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull", shape = 1))
        g2         g3 log(scale) 
    0.8466     1.2263     0.7579 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1916 0.1893 
> ## Cox
> sqrt(diag(vcov(cox)))
    g2     g3 
0.2089 0.2118 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1916 0.1893 
> 
> 
> 
> ## ************** Left-censored
> 
> mydata <- data.frame(y = Surv(y, sample(0:1, length(y), replace = TRUE), 
+                      type = "left"), g = g)
> 
> ## Estimate coefficients
> coef(opt <- mlt(m, data = mydata,  fixed = c("log(y)" = 1)))
     log(y) (Intercept)          g2          g3 
     1.0000      0.2164      0.9043      1.1407 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
    -0.2164     -0.9043     -1.1407 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1485 0.1474 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1916 0.1893 
> 
> 
> try(coef(cox <- coxph(y ~ g, data = mydata)))
Error in coxph(y ~ g, data = mydata) : 
  Cox model doesn't support "left" survival data
> try(coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull", shape = 1)))
Error in phreg(y ~ g, data = mydata, dist = "weibull", shape = 1) : 
  This model doesn't support "left" survival data
> 
> 
> 
> ## *************** Intervall-censored
> 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, fixed = c("log(y)" = 1)))
     log(y) (Intercept)          g2          g3 
     1.0000     -0.2613      0.7259      0.8381 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
     0.2613     -0.7259     -0.8381 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1651 0.1723 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1916 0.1893 
> 
> 
> try(coef(cox <- coxph(y ~ g, data = mydata)))
Error in coxph(y ~ g, data = mydata) : 
  Cox model doesn't support "interval" survival data
> try(coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull", shape = 1)))
Error in phreg(y ~ g, data = mydata, dist = "weibull", shape = 1) : 
  This model doesn't support "interval" survival data
> 
> 
> 
> 
> ## ************** Weibull - AFT *********************************
> 
> set.seed(196)
> 
> ### true dgp
> rY <- function(n, ...) rweibull(n, ...)
> pY <- function(x, ...) pweibull(x, ...)
> dY <- function(x, ...) dweibull(x, ...)
> 
> gf <- gl(3, 1)
> g <- rep(gf, 100)
> y <- rY(length(g), scale = (1:nlevels(g))[g], shape = 3)
> mydata <- data.frame(y = y, g = g)
> 
> boxplot(y ~ g, data = mydata)
> 
> Bb <- log_basis(numeric_var("y", support = range(y)), ui = "increasing", 
+                 remove_intercept = TRUE)
> Bx <- as.basis(~ g, data = mydata)
> m <- ctm(Bb, shifting = Bx, todist = "MinExtrVal")
> 
> ## Estimate coefficients
> 
> ## PH-scale
> (cf <- coef(opt <- mlt(m, data = mydata)))
     log(y) (Intercept)          g2          g3 
    3.10055    -0.02642    -2.12858    -3.37418 
> 
> (coef_cox <- coef(cox <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, 
+                  	          data = mydata)))
    g2     g3 
-2.114 -3.405 
> 
> (coef_phreg <- coef(phreg <- phreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g,
+                                    data = mydata, dist = "weibull")))
        g2         g3 log(scale) log(shape) 
  -2.12856   -3.37407    0.00851    1.13153 
> 
> ## AFT-scale
> coef(aft <- survreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata,
+                     dist = "weibull"))
(Intercept)          g2          g3 
    0.00851     0.68654     1.08827 
> 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
   -0.00852    -0.68652    -1.08825 
> coef_cox * aft$scale
     g2      g3 
-0.6817 -1.0984 
> coef_phreg[c("g2", "g3")] * aft$scale
     g2      g3 
-0.6865 -1.0883 
> 
> ## Compare shape parameters
> 1 / cf[1]
log(y) 
0.3225 
> 1 / exp(coef_phreg[c("log(shape)")])
log(shape) 
    0.3225 
> aft$scale
[1] 0.3225
> 
> ## Compare standard errors
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1695 0.2038 
> sqrt(diag(vcov(cox)))
    g2     g3 
0.1855 0.2196 
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1695 0.2038 
> 
> ## Compare log-Likelihoods
> logLik(aft)
'log Lik.' -259.1 (df=4)
> logLik(opt)
'log Lik.' -259.1 (df=4)
> 
> 
> ## ************************* Right-censored
> 
> mydata <- data.frame(y = Surv(y, sample(0:1, length(y), replace = TRUE)), g = g)
> 
> ## Estimate coefficients
> (cf <- coef(opt <- mlt(m, data = mydata)))
     log(y) (Intercept)          g2          g3 
     3.1005     -0.5711     -2.3824     -3.5226 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
    -0.1842     -0.7684     -1.1361 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
     0.1842      0.7684      1.1361 
> 
> coef_cox <- coef(cox <- coxph(y ~ g, data = mydata))
> coef_cox * aft$scale
     g2      g3 
-0.7957 -1.1705 
> 
> coefs_phreg <- coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull"))
> coefs_phreg[c("g2", "g3")] * aft$scale
     g2      g3 
-0.7684 -1.1361 
> 
> ## Compare standard errors
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.2378 0.2816 
> sqrt(diag(vcov(cox)))
    g2     g3 
0.2684 0.3115 
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.2378 0.2816 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -228.8 (df=4)
> logLik(aft)
'log Lik.' -228.8 (df=4)
> 
> 
> ## ****************** Left-censored
> 
> mydata <- data.frame(y = Surv(y, sample(0:1, length(y), replace = TRUE), 
+                      type = "left"), g = g)
> 
> ## Estimate coefficients
> (cf <- coef(opt <- mlt(m, data = mydata)))
     log(y) (Intercept)          g2          g3 
      2.190       0.583      -1.566      -2.389 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
     0.2662     -0.7151     -1.0906 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
    -0.2662      0.7151      1.0906 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -209.1 (df=4)
> logLik(aft)
'log Lik.' -209.1 (df=4)
> 
> 
> ## ************** Interval-censored
> 
> mydata <- data.frame(y = Surv(y, y + 1, sample(0:3, length(y), replace = TRUE),
+                      type = "interval"), g = g)
> 
> ## Estimate coefficients
> (cf <- coef(opt <- mlt(m, data = mydata)))
     log(y) (Intercept)          g2          g3 
     2.4686     -0.3285     -1.4644     -2.4999 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
    -0.1331     -0.5932     -1.0127 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
     0.1331      0.5932      1.0127 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -290.4 (df=4)
> logLik(aft)
'log Lik.' -290.4 (df=4)
> 
> 
> proc.time()
   user  system elapsed 
   2.51    0.29    2.79