Nothing
.frmt <- function(q) {
if (is.factor(q))
return(as.character(q))
return(formatC(q, digits = 3, width = 5))
}
### transformation function
tmlt <- function(object, newdata = NULL, q = NULL, log = FALSE, ...) {
if (log) warning("argument log ignored")
vn <- unlist(variable.names(object))
y <- variable.names(object, "response")
vnx <- vn[!(vn %in% y)]
model <- object$model
cf <- coef(object)
if ("bscaling" %in% names(model)) {
cf <- coef(object)
nd <- lapply(mkgrid(object, 2), function(x) x[1:2])
X <- model.matrix(model, data = as.data.frame(nd))
Assign <- attr(X, "Assign")
sterm <- exp(.5 * c(predict(model,
newdata = newdata,
coef = cf, terms = "bscaling")))
cf[Assign[2,] == "bscaling"] <- 0
cf <- matrix(cf, nrow = length(sterm), ncol = length(cf),
byrow = TRUE)
if (object$scale_shift) {
idx <- !Assign[2,] %in% "bscaling"
} else {
idx <- !Assign[2,] %in% c("bshifting", "bscaling")
}
cf[, idx] <- cf[, idx] * sterm
}
stopifnot(!is.null(newdata[[y]]) || !is.null(q))
if (is.data.frame(newdata)) {
### unconditional
if (length(vnx) == 0 & !is.null(q)) {
newdata <- data.frame(q)
names(newdata) <- y
q <- NULL
}
### in sample predictions
### this will _not_ work for censored responses
if (!is.null(newdata[[y]]) & is.null(q)) {
stopifnot(is.atomic(newdata[[y]]))
ret <- c(predict(model, newdata = newdata,
coef = cf, ...))
names(ret) <- rownames(newdata)
### P(Y \le y_K) = 1 but trafo can be < Inf
### depending on parameterisation
if (is.factor(f <- newdata[[y]])) {
i <- f == levels(f)[nlevels(f)]
if (any(i))
ret[i] <- Inf
} else {
### Y in (b[1], b[2]) => P(Y \le b[1]) = 0, P(Y \le b[2]) = 1
### <FIXME> what happens with deriv in ...? </FIXME>
b <- bounds(as.vars(object)[[y]])[[1]]
ret[f < (b[1] + .Machine$double.eps)] <- -Inf
ret[f > (b[2] - .Machine$double.eps)] <- Inf
}
return(ret)
}
### extra quantiles, compute transformation
### for each q and each row of newdata
stopifnot(is.atomic(q))
# stopifnot(length(unique(q)) == length(q))
dim <- c(length(q), nrow(newdata))
### <FIXME> this triggers a trick in
### basefun:::predict.basis; better checks needed </FIXME>
names(dim) <- c(y, vnx[1])
newdata <- as.list(newdata)
newdata[[y]] <- q
ret <- predict(model, newdata = newdata,
coef = cf, dim = dim, ...)
dn <- vector(mode = "list", length = 2)
names(dn) <- c(y, "newdata") ### deparse(substitute(newdata))) ?
dn[[y]] <- .frmt(q)
dn[["newdata"]] <- rownames(newdata)
dimnames(ret) <- dn
### trafo of last level is always Inf, see above
if (is.factor(f <- newdata[[y]])) {
i <- f == levels(f)[nlevels(f)]
if (any(i))
ret[i,] <- Inf
} else {
### Y in (b[1], b[2]) => P(Y \le b[1]) = 0, P(Y \le b[2]) = 1
### <FIXME> what happens with deriv in ...? </FIXME>
b <- bounds(as.vars(object)[[y]])[[1]]
ret[f < (b[1] + .Machine$double.eps)] <- -Inf
ret[f > (b[2] - .Machine$double.eps)] <- Inf
}
return(ret)
}
### need to generate newdata outside tmlt such that
### the rows of expand.grid(newdata) match the elements of
### the return value
stopifnot(is.atomic(newdata[[y]]))
stopifnot(is.null(q))
stopifnot(y %in% names(newdata))
ret <- predict(model, newdata = newdata,
coef = cf, dim = TRUE, ...)
dn <- lapply(newdata, .frmt)
dimnames(ret) <- dn
### trafo of last level is always Inf, see above
if (is.factor(f <- newdata[[y]])) {
i <- f == levels(f)[nlevels(f)]
if (any(i)) {
args <- lapply(names(dn), function(d) {
if (d == y)
return(i)
return(1:(dim(ret)[which(names(dn) == d)]))
})
ret <- do.call("[<-", c(list(i = ret), args,
list(value = Inf)))
}
} else {
### Y in (b[1], b[2]) => P(Y \le b[1]) = 0, P(Y \le b[2]) = 1
### <FIXME> what happens with deriv in ...? </FIXME>
b <- bounds(as.vars(object)[[y]])[[1]]
i <- (f < (b[1] + .Machine$double.eps))
if (any(i)) {
args <- lapply(names(dn), function(d) {
if (d == y)
return(i)
return(1:(dim(ret)[which(names(dn) == d)]))
})
ret <- do.call("[<-", c(list(i = ret), args,
list(value = -Inf)))
}
i <- (f > (b[2] - .Machine$double.eps))
if (any(i)) {
args <- lapply(names(dn), function(d) {
if (d == y)
return(i)
return(1:(dim(ret)[which(names(dn) == d)]))
})
ret <- do.call("[<-", c(list(i = ret), args,
list(value = Inf)))
}
}
return(ret)
}
### distribution function
pmlt <- function(object, newdata = NULL, q = NULL, lower.tail = TRUE, log = FALSE, ...)
object$todistr$p(tmlt(object = object, newdata = newdata,
q = q, ...), lower.tail = lower.tail, log.p = log)
### survivor function
smlt <- function(object, newdata = NULL, q = NULL, log = FALSE, ...)
pmlt(object = object, newdata = newdata, q = q, lower.tail = FALSE,
log = log, ...)
### cumulative hazard function
Hmlt <- function(object, newdata = NULL, q = NULL, log = FALSE, ...) {
if (object$todistr$name == "minimum extreme value") {
### Cox model
logH <- tmlt(object = object, newdata = newdata, q = q, ...)
if (log) return(logH)
return(exp(logH))
}
if (object$todistr$name == "exponential") {
### Aalen model
H <- tmlt(object = object, newdata = newdata, q = q, ...)
} else {
### generic
H <- -smlt(object = object, newdata = newdata, q = q, log = TRUE,
...)
}
if (log) return(log(H))
return(H)
}
### odds function
Omlt <- function(object, newdata = NULL, q = NULL, log = FALSE, ...) {
if (object$todistr$name == "logistic") {
logO <- tmlt(object = object, newdata = newdata, q = q, ...)
if (log) return(logO)
return(exp(logO))
}
F <- pmlt(object = object, newdata = newdata, q = q, log = log, ...)
S <- smlt(object = object, newdata = newdata, q = q, log = log, ...)
if (log) return(F - S)
return(F / S)
}
### numerically invert z = f(q)
.invf <- function(object, f, q, z) {
if (!is.matrix(f)) f <- matrix(f, nrow = 1)
N <- nrow(f)
K <- length(q)
stopifnot(ncol(f) == K)
if (!is.matrix(z)) z <- matrix(z, nrow = N)
stopifnot(nrow(z) == N)
nsim <- ncol(z)
y <- variable.names(object, "response")
bounds <- bounds(object)[[y]]
discrete <- !inherits(as.vars(object)[[y]],
"continuous_var")
if (discrete) {
### use "old" code
f <- f[rep(1:N, nsim), , drop = FALSE]
f <- cbind(-Inf, f, Inf)
i <- base::rowSums(f < as.vector(z))
return(q[i])
}
### use spline/approx to evaluate quantile function
nleft <- nexact <- nright <-
matrix(NA, nrow = N, ncol = nsim)
for (i in 1:N) {
pr <- f[i,]
pr0 <- which(pr < min(pr) + sqrt(.Machine$double.eps))
pr1 <- which(pr > max(pr) - sqrt(.Machine$double.eps))
rmi <- which(!is.finite(pr))
if (length(pr0) > 0L)
rmi <- c(rmi, pr0[-length(pr0)])
if (length(pr1) > 0L)
rmi <- c(rmi, pr1[-1L])
if (length(rmi) > 0L) {
qq <- q[-rmi]
pr <- pr[-rmi]
} else {
qq <- q
}
s <- spline(x = qq, y = pr, method = "hyman")
ynew <- approx(x = s$y, y = s$x, xout = z[i, ],
yleft = -Inf, yright = Inf)$y
if (bounds[1L] == min(q))
ynew[ynew == -Inf] <- min(q)
if (bounds[2L] == max(q))
ynew[ynew == Inf] <- max(q)
nleft[i, ynew == -Inf] <- bounds[1L]
nright[i, ynew == -Inf] <- min(q)
nleft[i, ynew == Inf] <- max(q)
nright[i, ynew == Inf] <- bounds[2L]
ynew[!is.finite(ynew)] <- NA
nexact[i, ] <- ynew
}
if (all(is.na(nleft)) && all(is.na(nright)))
return(as.vector(nexact))
return(R(as.vector(nexact),
cleft = as.vector(nleft),
cright = as.vector(nright)))
}
### quantile function
qmlt <- function(object, newdata = NULL, q = NULL, prob = .5, n = 50,
interpolate = FALSE, log = FALSE, ...) {
if (log) warning("argument log ignored")
if (interpolate)
warning("Argument interpolate ignored in mlt >= 1.2-1")
stopifnot(all(prob > sqrt(.Machine$double.eps)) &&
all(prob < 1 - sqrt(.Machine$double.eps)))
y <- variable.names(object, "response")
if (is.null(q))
q <- mkgrid(object, n = n)[[y]]
if (!is.null(newdata) & !is.data.frame(newdata)) {
newdata[[y]] <- NULL
nm <- names(newdata)
newdata[[y]] <- q
newdata <- newdata[c(y, nm)]
tr <- pmlt(object, newdata, ...)
} else {
tr <- pmlt(object, newdata = newdata, q = q, ...)
}
### convert potential array-valued distribution function
### to matrix where rows correspond to observations newdata
### and columns to quantiles q
trm <- matrix(tr, ncol = length(q), byrow = TRUE)
### Note: we invert the cdf direct;
### inverting transformation functions using approx didn't work well
qu <- matrix(prob,
nrow = nrow(trm), ncol = length(prob), byrow = TRUE)
ret <- .invf(object, f = trm, q = q, z = qu)
i <- as.vector(matrix(1:NROW(ret), ncol = nrow(trm), byrow = TRUE))
ret <- if (is.data.frame(ret)) ret[i,] else ret[i]
### arrays of factors are not allowed
if (is.factor(q)) return(ret)
### return "response" object
if (inherits(ret, "response"))
return(ret)
dim <- dim(tr)
dim[1] <- length(prob)
dn <- c(list(prob = .frmt(prob)), dimnames(tr)[-1L])
return(array(ret, dim = dim, dimnames = dn))
}
### density
dmlt <- function(object, newdata = NULL, q = NULL, log = FALSE, ...) {
y <- variable.names(object, "response")
response <- mkgrid(object, n = 10)[[y]]
### Lebesgue density only for double
if (.type_of_response(response) %in% c("double", "survival")) {
trafo <- tmlt(object, newdata = newdata, q = q, ...)
deriv <- 1
names(deriv) <- y
trafoprime <- tmlt(object, newdata = newdata, q = q,
deriv = deriv, ...)
### <FIXME> trafoprime is +/-Inf at boundaries, so use 0 density
trafoprime[!is.finite(trafoprime)] <- .Machine$double.eps
trafoprime <- pmax(.Machine$double.eps, trafoprime)
if (log)
return(object$todistr$d(trafo, log = TRUE) + log(trafoprime))
return(object$todistr$d(trafo) * trafoprime)
}
stopifnot(!is.null(newdata[[y]]) || !is.null(q))
### for factors and integers compute density as F(y) - F(y - 1)
lev <- levels(response)
if (is.data.frame(newdata)) {
### in sample density
if (!is.null(newdata[[y]]) & is.null(q)) {
stopifnot(is.atomic(newdata[[y]]))
q <- newdata[[y]]
first <- q == lev[1]
qwoK <- factor(lev[pmax(unclass(q) - 1, 1)],
levels = lev, labels = lev, ordered = is.ordered(q))
p <- pmlt(object, newdata = newdata, ...)
newdata[[y]] <- qwoK
pwoK <- pmlt(object, newdata = newdata, ...)
pwoK[first] <- 0
ret <- p - pwoK
} else {
### extra quantiles, compute density
### for each q and each row of newdata
stopifnot(is.atomic(q))
first <- q == lev[1]
qfirst <- q[first]
qwoK <- q[q != lev[length(lev)]]
qwo1 <- q[q != lev[1]]
pfirst <- pmlt(object, newdata = newdata, q = qfirst, ...)
pwo1 <- pmlt(object, newdata = newdata, q = qwo1, ...)
pwoK <- pmlt(object, newdata = newdata, q = qwoK, ...)
ret <- matrix(0, nrow = length(first), ncol = NCOL(pfirst))
ret[!first,] <- pwo1 - pwoK
ret[first,] <- pfirst
rownames(ret) <- as.character(q)
}
} else {
### need to generate newdata outside tmlt such that
### the rows of expand.grid(newdata) match the elements of
### the return value
stopifnot(is.atomic(newdata[[y]]))
stopifnot(is.null(q))
stopifnot(y %in% names(newdata))
dim <- sapply(newdata, NROW)
q <- newdata[[y]]
first <- q == lev[1]
qfirst <- q[first]
qwoK <- q[q != lev[length(lev)]]
qwo1 <- q[q != lev[1]]
newdata[[y]] <- qfirst
pfirst <- pmlt(object, newdata = newdata, ...)
newdata[[y]] <- qwo1
pwo1 <- pmlt(object, newdata = newdata, ...)
newdata[[y]] <- qwoK
pwoK <- pmlt(object, newdata = newdata, ...)
dn <- dim(pfirst)
names(dn) <- names(dimnames(pfirst))
frst <- lapply(names(dn), function(d) {
if (d != y) return(1:dn[d])
return(first)
})
ntfrst <- lapply(names(dn), function(d) {
if (d != y) return(1:dn[d])
return(!first)
})
dn <- lapply(names(dn), function(d) {
if (d != y) return(dimnames(pfirst)[[d]])
return(as.character(q))
})
ret <- array(0, dim = dim, dimnames = dn)
ret <- do.call("[<-",
c(list(i = ret), frst, list(value = pfirst)))
ret <- do.call("[<-",
c(list(i = ret), ntfrst, list(value = pwo1 - pwoK)))
}
if (log) return(log(ret))
return(ret)
}
### hazard function
hmlt <- function(object, newdata = object$data, q = NULL, log = FALSE, ...) {
y <- variable.names(object, "response")
response <- mkgrid(object, n = 10)[[y]]
if (.type_of_response(response) %in% c("double", "survival")) {
d <- dmlt(object, newdata = newdata, q = q, log = log, ...)
s <- smlt(object, newdata = newdata, q = q, log = log, ...)
if (log)
return(d - s)
return(d / s)
}
### discrete hazard: Prob(Y = y | Y >= y)
d <- dmlt(object, newdata = newdata, q = q, log = log, ...)
p <- pmlt(object, newdata = newdata, q = q, ...)
if (log)
return(d - log1p(-(p - exp(d))))
return(d / (1 - (p - d)))
}
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.