R Under development (unstable) (2023-12-17 r85691 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("tram") Loading required package: mlt Loading required package: basefun Loading required package: variables Loading required package: mvtnorm > library("multcomp") Loading required package: survival Loading required package: TH.data Loading required package: MASS Attaching package: 'TH.data' The following object is masked from 'package:MASS': geyser > > m <- BoxCox(dist ~ speed, data = cars) > m$negative [1] TRUE > coef(m) speed 0.3062335 > nd <- data.frame(speed = as.double(1:5 * 5)) > > PI(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.860529 0.9848214 0.9994191 0.9999926 1 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > PI(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.01517862 0.139471 0.5 0.860529 0.9848214 > PI(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.01517862 0.139471 0.5 0.860529 0.9848214 > PI(m, newdata = nd, reference = nd) 1 2 3 4 2 0.8605290 3 0.9848214 0.8605290 4 0.9994191 0.9848214 0.8605290 5 0.9999926 0.9994191 0.9848214 0.8605290 > > F <- m$model$todistr$p > Q <- m$model$todistr$q > f <- function(p, b) F(Q(p) - b) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.139471 with absolute error < 0.00012 > > OVL(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.443924 0.125728 0.02163296 0.002196175 0.0001292314 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > OVL(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.125728 0.443924 1 0.443924 0.125728 > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.125728 0.443924 1 0.443924 0.125728 > OVL(m, newdata = nd, reference = nd) 1 2 3 4 2 0.443924015 3 0.125727999 0.443924015 4 0.021632964 0.125727999 0.443924015 5 0.002196175 0.021632964 0.125727999 0.443924015 > > d <- m$model$todistr$d > Q <- m$model$todistr$q > f <- function(p, b) pmin(1, d(Q(p) - b) / d(Q(p))) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.4439239 with absolute error < 7.3e-05 > > m <- Colr(dist ~ speed, data = cars) > m$negative [1] FALSE > coef(m) speed -0.5293216 > nd <- data.frame(speed = as.double(1:5 * 5)) > > PI(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.8589555 0.9781801 0.9975257 0.9997579 0.9999781 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > PI(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.02181988 0.1410445 0.5 0.8589555 0.9781801 > PI(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.02181988 0.1410445 0.5 0.8589555 0.9781801 > PI(m, newdata = nd, reference = nd) 1 2 3 4 2 0.8589555 3 0.9781801 0.8589555 4 0.9975257 0.9781801 0.8589555 5 0.9997579 0.9975257 0.9781801 0.8589555 > > F <- m$model$todistr$p > Q <- m$model$todistr$q > f <- function(p, b) F(Q(p) + b) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.1410445 with absolute error < 9.2e-05 > > > OVL(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.4205382 0.1323967 0.03705084 0.01000088 0.002672582 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > OVL(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.1323967 0.4205382 1 0.4205382 0.1323967 > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.1323967 0.4205382 1 0.4205382 0.1323967 > OVL(m, newdata = nd, reference = nd) 1 2 3 4 2 0.42053821 3 0.13239674 0.42053821 4 0.03705084 0.13239674 0.42053821 5 0.01000088 0.03705084 0.13239674 0.42053821 > > d <- m$model$todistr$d > Q <- m$model$todistr$q > f <- function(p, b) pmin(1, d(Q(p) + b) / d(Q(p))) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.4205378 with absolute error < 3.4e-05 > > > m <- Coxph(dist ~ speed, data = cars) > m$negative [1] FALSE > coef(m) speed -0.280805 > nd <- data.frame(speed = as.double(1:5 * 5)) > > PI(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.8028218 0.9431093 0.9854007 0.9963744 0.9991071 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > PI(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.05689072 0.1971782 0.5 0.8028218 0.9431093 > PI(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.05689072 0.1971782 0.5 0.8028218 0.9431093 > PI(m, newdata = nd, reference = nd) 1 2 3 4 2 0.8028218 3 0.9431093 0.8028218 4 0.9854007 0.9431093 0.8028218 5 0.9963744 0.9854007 0.9431093 0.8028218 > > F <- m$model$todistr$p > Q <- m$model$todistr$q > f <- function(p, b) F(Q(p) + b) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.1971782 with absolute error < 2.2e-16 > > > OVL(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.5223833 0.2153214 0.07528464 0.02386655 0.007148041 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > OVL(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.2153214 0.5223833 1 0.5223833 0.2153214 > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.2153214 0.5223833 1 0.5223833 0.2153214 > OVL(m, newdata = nd, reference = nd) 1 2 3 4 2 0.52238327 3 0.21532141 0.52238327 4 0.07528464 0.21532141 0.52238327 5 0.02386655 0.07528464 0.21532141 0.52238327 > > d <- m$model$todistr$d > Q <- m$model$todistr$q > f <- function(p, b) pmin(1, d(Q(p) + b) / d(Q(p))) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.5223832 with absolute error < 2.8e-06 > > > m <- Lehmann(dist ~ speed, data = cars) > m$negative [1] TRUE > coef(m) speed 0.3478323 > nd <- data.frame(speed = as.double(1:5 * 5)) > > PI(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.8505805 0.9700647 0.9946083 0.9990486 0.9998327 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > PI(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.02993533 0.1494195 0.5 0.8505805 0.9700647 > PI(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.02993533 0.1494195 0.5 0.8505805 0.9700647 > PI(m, newdata = nd, reference = nd) 1 2 3 4 2 0.8505805 3 0.9700647 0.8505805 4 0.9946083 0.9700647 0.8505805 5 0.9990486 0.9946083 0.9700647 0.8505805 > > F <- m$model$todistr$p > Q <- m$model$todistr$q > f <- function(p, b) F(Q(p) - b) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.1494195 with absolute error < 6.7e-15 > > > OVL(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.4309584 0.1324664 0.03330627 0.007555076 0.001620911 > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > OVL(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.1324664 0.4309584 1 0.4309584 0.1324664 > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [,5] [1,] 0.1324664 0.4309584 1 0.4309584 0.1324664 > OVL(m, newdata = nd, reference = nd) 1 2 3 4 2 0.430958441 3 0.132466424 0.430958441 4 0.033306271 0.132466424 0.430958441 5 0.007555076 0.033306271 0.132466424 0.430958441 > > d <- m$model$todistr$d > Q <- m$model$todistr$q > f <- function(p, b) pmin(1, d(Q(p) - b) / d(Q(p))) > lp <- predict(m, newdata = nd, type = "lp") > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1]) 0.4309575 with absolute error < 5.2e-05 > > OVL(m, newdata = nd) [,1] [,2] [,3] [,4] [,5] [1,] 0.4309584 0.1324664 0.03330627 0.007555076 0.001620911 > OVL(m, newdata = nd, conf.level = .95) Estimate lwr upr 1 0.430958441 0.3206993334 0.56433049 2 0.132466424 0.0632140413 0.26158351 3 0.033306271 0.0098805486 0.10482119 4 0.007555076 0.0013940936 0.03827161 5 0.001620911 0.0001862683 0.01320327 attr(,"conf.level") [1] 0.95 attr(,"calpha") [1] 1.960371 > > lp15 <- c(predict(m, newdata = data.frame(speed = 15))) > OVL(m, newdata = nd, reference = lp15) [,1] [,2] [,3] [,4] [,5] [1,] 0.1324664 0.4309584 1 0.4309584 0.1324664 > OVL(m, newdata = nd, reference = lp15, conf.level = .95) Estimate lwr upr 1 0.1251517592 5.851127e-02 0.252039347 2 0.0312767073 9.088354e-03 0.100488045 3 0.0070724877 1.278423e-03 0.036576212 4 0.0015146035 1.705034e-04 0.012594431 5 0.0003135728 2.197845e-05 0.004187655 attr(,"conf.level") [1] 0.95 attr(,"calpha") [1] 1.960371 > > OVL(m, newdata = nd[-3,,drop = FALSE], reference = nd[3,,drop = FALSE]) [,1] [,2] [,3] [,4] [1,] 0.1324664 0.4309584 0.4309584 0.1324664 > OVL(m, newdata = nd[-3,,drop = FALSE], reference = nd[3,,drop = FALSE], conf.level = .95) Estimate lwr upr 3-1 0.1324664 0.06321515 0.2615796 3-2 0.4309584 0.32070160 0.5643272 3-4 0.4309584 0.32070160 0.5643272 3-5 0.1324664 0.06321515 0.2615796 attr(,"conf.level") [1] 0.95 attr(,"calpha") [1] 1.960326 > > OVL(m, newdata = nd, reference = nd) 1 2 3 4 2 0.430958441 3 0.132466424 0.430958441 4 0.033306271 0.132466424 0.430958441 5 0.007555076 0.033306271 0.132466424 0.430958441 > OVL(m, newdata = nd, reference = nd, conf.level = .95) Estimate lwr upr 1-2 0.430958441 0.320699820 0.56432978 1-3 0.132466424 0.063214279 0.26158268 2-3 0.430958441 0.320699820 0.56432978 1-4 0.033306271 0.009880609 0.10482062 2-4 0.132466424 0.063214279 0.26158268 3-4 0.430958441 0.320699820 0.56432978 1-5 0.007555076 0.001394105 0.03827131 2-5 0.033306271 0.009880609 0.10482062 3-5 0.132466424 0.063214279 0.26158268 4-5 0.430958441 0.320699820 0.56432978 attr(,"conf.level") [1] 0.95 attr(,"calpha") [1] 1.960362 > > OVL(m, newdata = nd, reference = nd, conf.level = .95, + calpha = univariate_calpha()) Estimate lwr upr 1-2 0.430958441 0.320719822 0.56430050 1-3 0.132466424 0.063224062 0.26154855 2-3 0.430958441 0.320719822 0.56430050 1-4 0.033306271 0.009883093 0.10479729 2-4 0.132466424 0.063224062 0.26154855 3-4 0.430958441 0.320719822 0.56430050 1-5 0.007555076 0.001394590 0.03825911 2-5 0.033306271 0.009883093 0.10479729 3-5 0.132466424 0.063224062 0.26154855 4-5 0.430958441 0.320719822 0.56430050 attr(,"conf.level") [1] 0.95 attr(,"calpha") [1] 1.959964 > OVL(m, newdata = nd, reference = nd, conf.level = .95, + calpha = adjusted_calpha()) Estimate lwr upr 1-2 0.430958441 0.320699820 0.56432978 1-3 0.132466424 0.063214279 0.26158268 2-3 0.430958441 0.320699820 0.56432978 1-4 0.033306271 0.009880609 0.10482062 2-4 0.132466424 0.063214279 0.26158268 3-4 0.430958441 0.320699820 0.56432978 1-5 0.007555076 0.001394105 0.03827131 2-5 0.033306271 0.009880609 0.10482062 3-5 0.132466424 0.063214279 0.26158268 4-5 0.430958441 0.320699820 0.56432978 attr(,"conf.level") [1] 0.95 attr(,"calpha") [1] 1.960362 > > proc.time() user system elapsed 2.54 0.32 2.82