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
> set.seed(29)
> tol <- sqrt(.Machine$double.eps)
> 
> n <- 100
> lcf <- c("(Intercept)" = .5, "g2" = 1, "x" = 2)
> sd <- .1
> cf <- c("y" = 1, -lcf) / sd
> cf <- cf[c("(Intercept)", "y", "g2", "x")]
> 
> d <- expand.grid(g = gl(2, 1), x = (1:n)/n)
> set.seed(29)
> d$y <- with(d, rnorm(nrow(d), 
+            mean = lcf["(Intercept)"] + lcf["g2"] * (g == "2") + lcf["x"] * x, 
+            sd = sd))
> 
> m <- ctm(polynomial_basis(numeric_var("y", support = range(d$y) * c(1, 1.1)),
+                             coef = c(TRUE, TRUE), ui = diag(2), ci = c(-Inf, 0)),
+            shift = ~ g + x, data = d)
> mod <- mlt(m, data = d, dofit = FALSE)
> coef(mod) <- cf
> 
> tfun <- function(d) with(expand.grid(d), 
+     y * cf["y"] + c(cf["(Intercept)"], sum(cf[c("(Intercept)", "g2")]))[g] + cf["x"] * x)
> pfun <- function(d) with(expand.grid(d), pnorm(y, 
+            mean = lcf["(Intercept)"] + lcf["g2"] * (g == "2") + lcf["x"] * x,
+            sd = sd))
> dfun <- function(d) with(expand.grid(d), dnorm(y, 
+            mean = lcf["(Intercept)"] + lcf["g2"] * (g == "2") + lcf["x"] * x,
+            sd = sd))
> qfun <- function(d) with(expand.grid(d), qnorm(p,
+            mean = lcf["(Intercept)"] + lcf["g2"] * (g == "2") + lcf["x"] * x,
+            sd = sd))
> 
> (ny <- mkgrid(mod, 10)$y)
 [1] 0.3916629 0.7825718 1.1734807 1.5643897 1.9552986 2.3462075 2.7371164
 [8] 3.1280254 3.5189343 3.9098432
> nd <- list(y = ny, g = gl(2, 1), x = (1:10) / 10)
> ndx <- expand.grid(nd[-1])
> end <- ndx
> end$y <- seq(from = min(d$y), to = max(d$y), length = nrow(end))
> 
> max(abs(predict(mod, newdata = end) - model.matrix(m, data = end) %*% cf)) < tol
[1] TRUE
> 
> cf2 <- cf
> cf2[1:2] <- 0
> max(abs(predict(mod, newdata = end, terms = "bshifting") - model.matrix(m, data = end) %*% cf2)) < tol
[1] TRUE
> 
> p1 <- predict(mod, newdata = ndx, q = ny)
> 
> p2 <- predict(mod, newdata = nd)
> 
> max(abs(c(p1) - c(p2))) < tol
[1] TRUE
> 
> p3 <- tfun(nd)
> max(abs(c(p2) - c(p3))) < tol
[1] TRUE
> 
> p1 <- predict(mod, newdata = nd, type = "distribution")
> p2 <- pfun(nd)
> max(abs(c(p1) - c(p2))) < tol
[1] TRUE
> 
> p1 <- predict(mod, newdata = nd, type = "density")
> p2 <- dfun(nd)
> max(abs(c(p1) - c(p2))) < tol
[1] TRUE
> 
> p1 <- predict(mod, newdata = nd, type = "quantile", prob = 1:9 / 10, K = 50000)
> 
> min(apply(matrix(p1, nrow = 9), 2, diff))
[1] 0.02533471
> 
> ### more than one stratifying variable
> n <- 5
> nd <- expand.grid(d <- list(s1 = gl(3, n), s2 = gl(3, n), x = runif(n)))
> nd$y <- rnorm(nrow(nd), mean = (1:3)[nd$s1] + 2 * nd$x, sd = (1:3)[nd$s2]/3)
> yvar <- numeric_var("y", support = quantile(nd$y, prob = c(2, 8)/10))
> b1 <- as.basis(~ s1 - 1, data = nd)
> b2 <- as.basis(~ s2 - 1, data = nd)
> mctm <- ctm(Bernstein_basis(yvar, order = 2, ui = "increasing"),
+             interacting = b(b1 = b1, b2 = b2),
+             shifting = ~x, data = nd)
> m <- mlt(mctm, data = nd, scale = TRUE)
> ## in-sample predictions (for given response)
> p1 <- predict(m)
> p2 <- predict(m, newdata = nd)
> stopifnot(isTRUE(all.equal(max(abs(c(p1) - c(p2))), 0)))
> ## evaluate model for grid of response values
> pnd <- nd
> pnd$y <- NULL
> p1 <- predict(m, newdata = pnd, K = 21)
> d$y <- mkgrid(m, n = 21)[["y"]]
> d <- d[c("y", "s1", "s2", "x")]
> p2 <- predict(m, newdata = d)
> stopifnot(isTRUE(all.equal(max(abs(c(p1) - c(p2))), 0)))
> 
> type <- c("distribution", "density", "survivor", "hazard", "cumhazard", "odds")
> 
> out <- sapply(type, function(ty) {
+  p1 <- log(predict(m, newdata = pnd, type = ty))
+  p2 <- predict(m, newdata = pnd, type = ty, log = TRUE)
+  p3 <- predict(m, newdata = pnd, type = paste0("log", ty))
+  stopifnot(isTRUE(all.equal(p1, p2)))
+  stopifnot(isTRUE(all.equal(p1, p3)))
+ })
> 
> 
> proc.time()
   user  system elapsed 
   3.70    0.75    4.43