Nothing
library("mlt")
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)
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
cf2 <- cf
cf2[1:2] <- 0
max(abs(predict(mod, newdata = end, terms = "bshifting") - model.matrix(m, data = end) %*% cf2)) < tol
p1 <- predict(mod, newdata = ndx, q = ny)
p2 <- predict(mod, newdata = nd)
max(abs(c(p1) - c(p2))) < tol
p3 <- tfun(nd)
max(abs(c(p2) - c(p3))) < tol
p1 <- predict(mod, newdata = nd, type = "distribution")
p2 <- pfun(nd)
max(abs(c(p1) - c(p2))) < tol
p1 <- predict(mod, newdata = nd, type = "density")
p2 <- dfun(nd)
max(abs(c(p1) - c(p2))) < tol
p1 <- predict(mod, newdata = nd, type = "quantile", prob = 1:9 / 10, K = 50000)
min(apply(matrix(p1, nrow = 9), 2, diff))
### 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)))
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.