R Under development (unstable) (2025-11-17 r89032 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.1190 
> 
> 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.07598     0.72341     1.02746 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
   -0.08224     0.78304     1.11216 
> 
> 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.6851      0.7375      1.1013 
> 
> ## Estimate coefficients 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
     0.6851     -0.7374     -1.1014 
> coef(cox <- coxph(y ~ g, data = mydata))
    g2     g3 
0.5922 0.9482 
> coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull", shape = 1))
        g2         g3 log(scale) 
    0.7374     1.1014     0.6851 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1899 0.1882 
> ## Cox
> sqrt(diag(vcov(cox)))
    g2     g3 
0.2024 0.2052 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1899 0.1882 
> 
> 
> 
> ## ************** 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.2317      1.0233      1.0506 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
    -0.2317     -1.0233     -1.0506 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1510 0.1466 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1899 0.1882 
> 
> 
> 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.3657      0.6942      1.1192 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
     0.3657     -0.6942     -1.1192 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.1706 0.1719 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.1899 0.1882 
> 
> 
> 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.10042    -0.02639    -2.12857    -3.37410 
> 
> (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.008512   -0.686543   -1.088271 
> 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.2322     -0.7415     -2.0683     -3.5958 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
    -0.2294     -0.6399     -1.1125 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
     0.2294      0.6399      1.1125 
> 
> coef_cox <- coef(cox <- coxph(y ~ g, data = mydata))
> coef_cox * aft$scale
     g2      g3 
-0.6306 -1.0982 
> 
> coefs_phreg <- coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull"))
> coefs_phreg[c("g2", "g3")] * aft$scale
     g2      g3 
-0.6399 -1.1125 
> 
> ## Compare standard errors
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
    g2     g3 
0.2349 0.2931 
> sqrt(diag(vcov(cox)))
    g2     g3 
0.2517 0.3065 
> sqrt(diag(phreg$var))[c("g2", "g3")]
    g2     g3 
0.2349 0.2931 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -226.3 (df=4)
> logLik(aft)
'log Lik.' -226.3 (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.1517      0.4396     -1.2440     -2.2483 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
     0.2043     -0.5782     -1.0449 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
    -0.2043      0.5782      1.0449 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -227 (df=4)
> logLik(aft)
'log Lik.' -227 (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.474      -0.366      -1.466      -2.405 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
    -0.1479     -0.5927     -0.9723 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
     0.1479      0.5927      0.9722 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -291.4 (df=4)
> logLik(aft)
'log Lik.' -291.4 (df=4)
> 
> 
> proc.time()
   user  system elapsed 
   2.48    0.09    2.56