R Under development (unstable) (2023-06-29 r84618 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("mlt") Loading required package: basefun Loading required package: variables > library("numDeriv") > set.seed(29) > options(digits = 4) > > n <- 20 > ### just for interface checking > ### we need something better! > d <- data.frame(x1 = 1:n, x2 = sample(1:n) + 1, y = rnorm(n)) > m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y)), + coef = c(TRUE, TRUE), ci = c(-Inf, 0)), + shift = ~ x1 + x2, data = d) > mod <- mlt(m, data = d) > coef(mod) (Intercept) y x1 x2 -1.15423 1.11402 0.02335 0.06000 > > p <- predict(mod, newdata = d) > > K <- 15 > q <- mkgrid(m, n = K)[["y"]] > p1 <- predict(mod, newdata = d[, c("x1", "x2")], q = q) > p2 <- predict(mod, newdata = d[, c("x1", "x2")], K = K) > stopifnot(all.equal(p1, p2)) > > p0 <- predict(mod$model$model, + newdata = expand.grid(d), coef = coef(mod)) > p1 <- predict(mod, newdata = as.list(d)) > p2 <- predict(mod, newdata = d, q = d$y[1]) > > max(abs(p0 - as.vector(p1))) [1] 0 > > all.equal(p1[cbind(1:n, 1:n, 1), drop = TRUE], + drop(p2)) [1] TRUE > > all.equal(p1[cbind(1:n, 1:n, 1:n), drop = TRUE], + drop(p), check.attributes = FALSE) [1] TRUE > > predict(mod, newdata = list(x1 = 1:3, x2 = 2:3), prob = c(.25, .5), type = "quantile") , , x2 = 2 x1 prob 1 2 3 0.25 0.3019 0.2810 0.2600 0.5 0.9074 0.8865 0.8655 , , x2 = 3 x1 prob 1 2 3 0.25 0.2481 0.2271 0.2061 0.5 0.8536 0.8326 0.8116 > > simulate(mod, nsim = 1, seed = 291, interpolate = FALSE) [1] -0.98947 0.97474 -0.78103 1.58963 [5] 0.16146 1.04445 (2.198, Inf] 1.56490 [9] -0.37026 0.85677 -0.59923 -0.30236 [13] 0.65478 0.34240 1.01619 1.06508 [17] ( -Inf, -1.449] 0.07819 -0.17235 1.46074 > > d$y <- gl(3, 1, ordered = TRUE)[rep(1:3, length = n)] > > r <- as.basis(d$y) #as.basis(~ y, data = d, remove_intercept = TRUE, > # contrasts.arg = list(y = function(n) > # contr.treatment(n, base = 3)), > # ui = diff(diag(2)), ci = 0) > > mod2 <- mlt(ctm(r, shift = ~ x1 + x2, data = d), data = d) > > predict(mod2, q = unique(d$y)) y [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] 1 -0.1534 -0.4301 -0.4719 -0.04421 -0.7513 -0.08881 -0.2481 -0.2899 -0.7231 2 0.7796 0.5030 0.4611 0.88882 0.1817 0.84423 0.6850 0.6431 0.2099 3 Inf Inf Inf Inf Inf Inf Inf Inf Inf y [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] 1 -0.6476 -0.0242 -0.6922 -0.4210 -0.3454 -0.1525 -0.2335 -0.6276 -0.552 2 0.2854 0.9088 0.2409 0.5121 0.5876 0.7805 0.6995 0.3055 0.381 3 Inf Inf Inf Inf Inf Inf Inf Inf Inf y [,19] [,20] 1 -0.5939 -0.4010 2 0.3391 0.5321 3 Inf Inf > > predict(mod2, prob = 1:9 / 10, type = "quantile") [1] 1 1 1 1 2 2 2 3 3 1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3 1 1 1 1 2 2 2 2 3 1 [38] 1 2 2 2 3 3 3 3 1 1 1 1 2 2 2 2 3 1 1 1 1 2 2 2 3 3 1 1 1 2 2 2 2 3 3 1 1 [75] 2 2 2 3 3 3 3 1 1 2 2 2 2 3 3 3 1 1 1 1 2 2 2 2 3 1 1 2 2 2 3 3 3 3 1 1 1 [112] 2 2 2 3 3 3 1 1 1 2 2 2 2 3 3 1 1 1 1 2 2 2 3 3 1 1 1 1 2 2 2 3 3 1 1 2 2 [149] 2 2 3 3 3 1 1 2 2 2 2 3 3 3 1 1 2 2 2 2 3 3 3 1 1 1 2 2 2 2 3 3 Levels: 1 < 2 < 3 > > simulate(mod2, nsim = 3, seed = 29) [[1]] [1] 1 1 1 1 3 1 3 3 1 1 3 2 1 2 1 3 3 2 3 2 Levels: 1 < 2 < 3 [[2]] [1] 2 2 3 3 3 1 2 3 3 1 1 1 1 3 2 1 3 2 2 1 Levels: 1 < 2 < 3 [[3]] [1] 2 1 3 2 1 3 3 2 3 3 1 3 1 2 2 2 2 3 2 1 Levels: 1 < 2 < 3 > > predict(mod2, q = unique(d$y), type = "density") [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] 1 0.4390 0.3336 0.3185 0.4824 0.2262 0.4646 0.4020 0.3859 0.2348 0.2586 0.4903 2 0.3432 0.3589 0.3592 0.3306 0.3459 0.3361 0.3513 0.3540 0.3483 0.3537 0.3279 3 0.2178 0.3075 0.3224 0.1870 0.4279 0.1993 0.2467 0.2601 0.4169 0.3877 0.1817 [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] 1 0.2444 0.3369 0.3649 0.4394 0.4077 0.2651 0.2905 0.2763 0.3442 2 0.3508 0.3588 0.3567 0.3431 0.3502 0.3549 0.3579 0.3565 0.3584 3 0.4048 0.3043 0.2784 0.2175 0.2421 0.3800 0.3516 0.3673 0.2973 > > predict(mod2, list(y = unique(d$y), x1 = 1:3, x2 = 2:3), type = "density") , , 2 1 2 3 1 0.5012 0.5002 0.4991 2 0.3242 0.3245 0.3249 3 0.1746 0.1753 0.1760 , , 3 1 2 3 1 0.4856 0.4845 0.4835 2 0.3295 0.3299 0.3302 3 0.1849 0.1856 0.1863 > > ### some basis checks: continuous > d <- data.frame(x1 = 1:n, x2 = sample(1:n) + 1, y = rnorm(n)) > m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y)), + coef = c(TRUE, TRUE), ci = c(-Inf, 0)), + shift = ~ x1 + x2, data = d) > mod <- mlt(m, data = d) > > .chk <- function(x) + stopifnot(isTRUE(max(abs(x), na.rm = TRUE) < sqrt(.Machine$double.eps))) > > cont <- quote({ + nd <- d + nd$y <- NULL + q <- mkgrid(mod, 10)[[1]] + p <- predict(mod, newdata = nd, q = q, type = "distribution") + s <- predict(mod, newdata = nd, q = q, type = "survivor") + .chk(predict(mod, newdata = nd, q = q, type = "distribution", log = TRUE) - log(p)) + .chk(predict(mod, newdata = nd, q = q, type = "distribution", lower.tail = FALSE) - s) + .chk(predict(mod, newdata = nd, q = q, type = "distribution", lower.tail = + FALSE, log = TRUE) - log(s)) + + o <- predict(mod, newdata = nd, q = q, type = "odds") + .chk(o - p / s) + + df <- function(q) + predict(mod, newdata = nd[1,], q = q, type = "distribution") + da <- sapply(q, function(q) grad(df, q)) + dd <- predict(mod, newdata = nd[1,], q = q, type = "density") + .chk(da - dd) + + h <- predict(mod, newdata = nd, q = q, type = "hazard") + .chk(predict(mod, newdata = nd, q = q, type = "loghazard") - log(h)) + + H <- predict(mod, newdata = nd, q = q, type = "cumhazard") + .chk(H + log(s)) + + .chk(predict(mod, newdata = nd, q = q, type = "logcumhazard") - log(H)) + + dh <- function(q) + predict(mod, newdata = nd[1,], q = q, type = "cumhazard") + + da <- sapply(q, function(q) grad(dh, q)) + dd <- predict(mod, newdata = nd[1,], q = q, type = "hazard") + .chk(da - dd) + }) > > m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y)), + coef = c(TRUE, TRUE), ci = c(-Inf, 0)), + shift = ~ x1 + x2, data = d, todistr = "Normal") > mod <- mlt(m, data = d) > > eval(cont) > > m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y)), + coef = c(TRUE, TRUE), ci = c(-Inf, 0)), + shift = ~ x1 + x2, data = d, todistr = "Logistic") > mod <- mlt(m, data = d) > > eval(cont) > > m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y)), + coef = c(TRUE, TRUE), ci = c(-Inf, 0)), + shift = ~ x1 + x2, data = d, todistr = "MinExtrVal") > mod <- mlt(m, data = d) > > eval(cont) > > m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y)), + coef = c(TRUE, TRUE), ci = c(-Inf, 0)), + shift = ~ x1 + x2, data = d, todistr = "MaxExtrVal") > mod <- mlt(m, data = d) > > eval(cont) > > ### some basis checks: discrete > > disc <- quote({ + nd <- d + nd$y <- NULL + q <- mkgrid(mod, 10)[[1]] + p <- predict(mod, newdata = nd, q = q, type = "distribution") + s <- predict(mod, newdata = nd, q = q, type = "survivor") + .chk(predict(mod, newdata = nd, q = q, type = "distribution", log = TRUE) - log(p)) + .chk(predict(mod, newdata = nd, q = q, type = "distribution", lower.tail = FALSE) - s) + .chk(predict(mod, newdata = nd, q = q, type = "distribution", lower.tail = + FALSE, log = TRUE) - log(s)) + + o <- predict(mod, newdata = nd, q = q, type = "odds") + .chk(o - p / s) + + dd <- predict(mod, newdata = nd, q = q, type = "density") + + .chk(apply(p, 2, function(x) diff(c(0, x))) - dd) + + h <- predict(mod, newdata = nd, q = q, type = "hazard") + + .chk(dd / (1 - (p - dd)) - h) + + .chk(apply(h, 2, function(x) cumprod(1 - x)) - s) + + H <- predict(mod, newdata = nd, q = q, type = "cumhazard") + + .chk(H + log(s)) + }) > > d <- data.frame(x1 = 1:n, x2 = sample(1:n) + 1, y = sample(gl(4, 5, ordered = TRUE))) > m <- ctm(as.basis(d$y), shift = ~ x1 + x2, data = d, todistr = "Normal") > mod <- mlt(m, data = d) > > eval(disc) > > d <- data.frame(x1 = 1:n, x2 = sample(1:n) + 1, y = sample(gl(4, 5, ordered = TRUE))) > m <- ctm(as.basis(d$y), shift = ~ x1 + x2, data = d, todistr = "Logistic") > mod <- mlt(m, data = d) > > eval(disc) > > d <- data.frame(x1 = 1:n, x2 = sample(1:n) + 1, y = sample(gl(4, 5, ordered = TRUE))) > m <- ctm(as.basis(d$y), shift = ~ x1 + x2, data = d, todistr = "MinExtrVal") > mod <- mlt(m, data = d) > > eval(disc) > > d <- data.frame(x1 = 1:n, x2 = sample(1:n) + 1, y = sample(gl(4, 5, ordered = TRUE))) > m <- ctm(as.basis(d$y), shift = ~ x1 + x2, data = d, todistr = "MaxExtrVal") > mod <- mlt(m, data = d) > > eval(disc) > > > proc.time() user system elapsed 7.59 0.50 8.06