R Under development (unstable) (2024-11-18 r87347 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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)
> 
> ## ************** 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.0000000  -0.1230066   0.7915161   1.1190543 
> 
> coef(aft <- survreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata,
+                     dist = "exponential"))
(Intercept)          g2          g3 
  0.1230066  -0.7915161  -1.1190543 
> 
> coef(cox <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata))
      g2       g3 
0.699808 1.022699 
> 
> coef(phreg <- phreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata, 
+                     dist = "weibull", shape = 1))
        g2         g3 log(scale) 
 0.7915161  1.1190543  0.1230066 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
       g2        g3 
0.1414214 0.1414214 
> ## Cox
> sqrt(diag(vcov(cox)))
       g2        g3 
0.1523056 0.1557383 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
       g2        g3 
0.1414214 0.1414214 
> 
> ## Compare log-Likelihoods
> logLik(aft)
'log Lik.' -145.845 (df=3)
> logLik(opt)
'log Lik.' -145.845 (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.92384441 -0.07596598  0.72339973  1.02745353 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
-0.08222811  0.78303199  1.11214996 
> 
> coef(aft2 <- survreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata,
+                      dist = "weibull"))
(Intercept)          g2          g3 
  0.0822325  -0.7830360  -1.1121547 
> 
> ## Compare Weibull shape paramters
> 1 / cf[1]
  log(y) 
1.082433 
> aft2$scale
[1] 1.082433
> 
> ## Compare log-Likelihoods
> logLik(opt2)
'log Lik.' -144.2832 (df=4)
> logLik(aft2)
'log Lik.' -144.2832 (df=4)
> 
> sqrt(diag(vcov(opt2)))[c("g2", "g3")]
       g2        g3 
0.1462995 0.1502459 
> sqrt(diag(vcov(aft2)))[c("g2", "g3")]
       g2        g3 
0.1531531 0.1531334 
> 
> 
> ## *************** 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.0000000  -0.7578848   0.8465755   1.2262994 
> 
> ## Estimate coefficients 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
  0.7578849  -0.8465759  -1.2262999 
> coef(cox <- coxph(y ~ g, data = mydata))
       g2        g3 
0.8332798 1.2134218 
> coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull", shape = 1))
        g2         g3 log(scale) 
 0.8465759  1.2262999  0.7578849 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
       g2        g3 
0.1916379 0.1892540 
> ## Cox
> sqrt(diag(vcov(cox)))
       g2        g3 
0.2088555 0.2118118 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
       g2        g3 
0.1916379 0.1892540 
> 
> 
> 
> ## ************** 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.0000000   0.2164302   0.9043383   1.1406712 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
 -0.2164302  -0.9043383  -1.1406712 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
       g2        g3 
0.1485446 0.1473571 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
       g2        g3 
0.1916379 0.1892540 
> 
> 
> 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.0000000  -0.2612821   0.7258507   0.8381451 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "expo"))
(Intercept)          g2          g3 
  0.2612821  -0.7258509  -0.8381449 
> 
> ## Compare standard errors
> ## MLT
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
       g2        g3 
0.1651103 0.1723239 
> ## phreg
> sqrt(diag(phreg$var))[c("g2", "g3")]
       g2        g3 
0.1916379 0.1892540 
> 
> 
> 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.10055012 -0.02641571 -2.12858193 -3.37418384 
> 
> (coef_cox <- coef(cox <- coxph(Surv(y, rep(TRUE, nrow(mydata))) ~ g, 
+                  	          data = mydata)))
       g2        g3 
-2.113511 -3.405349 
> 
> (coef_phreg <- coef(phreg <- phreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g,
+                                    data = mydata, dist = "weibull")))
         g2          g3  log(scale)  log(shape) 
-2.12855989 -3.37407494  0.00851049  1.13153109 
> 
> ## AFT-scale
> coef(aft <- survreg(Surv(y, rep(TRUE, nrow(mydata))) ~ g, data = mydata,
+                     dist = "weibull"))
(Intercept)          g2          g3 
 0.00851049  0.68654367  1.08827090 
> 
> cf[-1] / cf[1]
 (Intercept)           g2           g3 
-0.008519685 -0.686517504 -1.088253279 
> coef_cox * aft$scale
        g2         g3 
-0.6816897 -1.0983578 
> coef_phreg[c("g2", "g3")] * aft$scale
        g2         g3 
-0.6865437 -1.0882709 
> 
> ## Compare shape parameters
> 1 / cf[1]
   log(y) 
0.3225234 
> 1 / exp(coef_phreg[c("log(shape)")])
log(shape) 
  0.322539 
> aft$scale
[1] 0.322539
> 
> ## Compare standard errors
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
       g2        g3 
0.1695244 0.2038101 
> sqrt(diag(vcov(cox)))
       g2        g3 
0.1855397 0.2196354 
> sqrt(diag(phreg$var))[c("g2", "g3")]
       g2        g3 
0.1695251 0.2038084 
> 
> ## Compare log-Likelihoods
> logLik(aft)
'log Lik.' -259.0559 (df=4)
> logLik(opt)
'log Lik.' -259.0559 (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.1004541  -0.5711125  -2.3823821  -3.5225652 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
 -0.1842029  -0.7683978  -1.1361449 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
  0.1842060   0.7683967   1.1361412 
> 
> coef_cox <- coef(cox <- coxph(y ~ g, data = mydata))
> coef_cox * aft$scale
        g2         g3 
-0.7957308 -1.1705361 
> 
> coefs_phreg <- coef(phreg <- phreg(y ~ g, data = mydata, dist = "weibull"))
> coefs_phreg[c("g2", "g3")] * aft$scale
        g2         g3 
-0.7683967 -1.1361412 
> 
> ## Compare standard errors
> sqrt(diag(vcov(opt)))[c("g2", "g3")]
       g2        g3 
0.2377624 0.2816145 
> sqrt(diag(vcov(cox)))
       g2        g3 
0.2684431 0.3114552 
> sqrt(diag(phreg$var))[c("g2", "g3")]
       g2        g3 
0.2377627 0.2816141 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -228.7986 (df=4)
> logLik(aft)
'log Lik.' -228.7986 (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.1902538   0.5830063  -1.5662463  -2.3886791 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
  0.2661821  -0.7150981  -1.0905947 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
 -0.2661692   0.7150879   1.0905913 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -209.1118 (df=4)
> logLik(aft)
'log Lik.' -209.1118 (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.4686323  -0.3284947  -1.4644029  -2.4999036 
> cf[-1] / cf[1]
(Intercept)          g2          g3 
 -0.1330675  -0.5932041  -1.0126675 
> 
> coef(aft <- survreg(y ~ g, data = mydata, dist = "weibull"))
(Intercept)          g2          g3 
  0.1330675   0.5932042   1.0126677 
> 
> ## Compare log-Likelihoods
> logLik(opt)
'log Lik.' -290.4319 (df=4)
> logLik(aft)
'log Lik.' -290.4319 (df=4)
> 
> 
> proc.time()
   user  system elapsed 
   3.59    0.37    3.95