Nothing
.drop <- function(x) {
if (inherits(x, "Matrix")) x <- x@x
drop(x)
}
.pmax <- function(x, y) {
y[y < x] <- x
return(y)
}
.xmb <- function(X, beta) {
if (is.matrix(beta)) {
stopifnot(NROW(beta) == NROW(X))
return(.Call("R_rowSums", X, beta))
# return(base::rowSums(X * beta))
}
return(X %*% beta)
}
### if (finite) fun(X %*% beta + offset) * Xmult else value
.dealinf <- function(X, beta, offset, fun, value, Xmult = FALSE) {
if (is.null(X)) return(value)
lpr <- .drop(offset + .xmb(X, beta))
OK <- is.finite(lpr)
if (all(OK)) {
ret <- fun(lpr)
} else {
ret <- rep(value, length(lpr))
if (all(!OK)) {
if (!Xmult) return(ret)
return(matrix(ret, nrow = nrow(X), ncol = ncol(X)))
}
ret[OK] <- fun(lpr[OK])
}
if (Xmult) {
### this is costly, can we avoid this?
if (!all(OK))
X[!OK,] <- 0
ret <- ret * X
}
return(ret)
}
.log <- function(x)
return(log(.pmax(.Machine$double.eps, x)))
..mlt_loglik_interval <- function(d, mml, mmr, offset = 0, beta) {
RC <- !is.finite(mmr[,1])
LC <- !is.finite(mml[,1])
ret <- numeric(nrow(mmr))
### right-censored: ret = log(1 - F(...))
if (any(irc <- RC & !LC)) {
if (is.matrix(beta))
ret[irc] <- .dealinf(mml[irc,,drop = FALSE], beta[irc,,drop = FALSE],
offset[irc], function(x) d$p(x, lower.tail = FALSE, log.p = TRUE), 0)
else
ret[irc] <- .dealinf(mml[irc,,drop = FALSE], beta, offset[irc],
function(x) d$p(x, lower.tail = FALSE, log.p = TRUE), 0)
}
### left-censored: ret = log(F(...))
if (any(ilc <- !RC & LC)) {
if (is.matrix(beta))
ret[ilc] <- .dealinf(mmr[ilc,,drop = FALSE], beta[ilc,,drop = FALSE],
offset[ilc], function(x) d$p(x, lower.tail = TRUE, log.p = TRUE), 1)
else
ret[ilc] <- .dealinf(mmr[ilc,,drop = FALSE], beta, offset[ilc],
function(x) d$p(x, lower.tail = TRUE, log.p = TRUE), 1)
}
### interval-censored: log(F(...) - F(...))
iic <- !RC & !LC
if (is.matrix(beta))
ret[iic] <- .log(.dealinf(mmr[iic,,drop = FALSE], beta[iic,,drop = FALSE],
offset[iic], d$p, 1) -
.dealinf(mml[iic,,drop = FALSE], beta[iic,,drop = FALSE],
offset[iic], d$p, 0))
else
ret[iic] <- .log(.dealinf(mmr[iic,,drop = FALSE], beta, offset[iic], d$p, 1) -
.dealinf(mml[iic,,drop = FALSE], beta, offset[iic], d$p, 0))
ret
}
.mlt_loglik_interval <- function(d, mml, mmr, offset = 0, mmtrunc = NULL) {
if (is.null(mmtrunc))
return(function(beta)
..mlt_loglik_interval(d, mml, mmr, offset, beta))
return(function(beta)
..mlt_loglik_interval(d, mml, mmr, offset, beta) -
..mlt_loglik_interval(d, mmtrunc$left,
mmtrunc$right, offset, beta))
}
..mlt_score_interval <- function(d, mml, mmr, offset = 0, beta, Xmult = TRUE)
(.dealinf(mmr, beta, offset, d$d, 0, Xmult = Xmult) -
.dealinf(mml, beta, offset, d$d, 0, Xmult = Xmult)) /
.pmax(.Machine$double.eps^(1/3), (.dealinf(mmr, beta, offset, d$p, 1) -
.dealinf(mml, beta, offset, d$p, 0)))
.mlt_score_interval <- function(d, mml, mmr, offset = 0, mmtrunc = NULL) {
if (is.null(mmtrunc))
return(function(beta, Xmult = TRUE)
..mlt_score_interval(d, mml, mmr, offset, beta, Xmult))
return(function(beta, Xmult = TRUE)
..mlt_score_interval(d, mml, mmr, offset, beta, Xmult) -
..mlt_score_interval(d, mmtrunc$left,
mmtrunc$right, offset, beta, Xmult))
}
..mlt_hessian_interval <- function(d, mml, mmr, offset = 0, w = 1, beta) {
Fr <- .dealinf(mmr, beta, offset, d$p, 1)
Fl <- .dealinf(mml, beta, offset, d$p, 0)
fr <- .dealinf(mmr, beta, offset, d$d, 0)
fl <- .dealinf(mml, beta, offset, d$d, 0)
dfr <- .dealinf(mmr, beta, offset, d$dd, 0)
dfl <- .dealinf(mml, beta, offset, d$dd, 0)
if (length(w) != length(Fr)) w <- rep(w, length(Fr))
Frl <- Fr - Fl
w1 <- dfr / Frl * w
w2 <- dfl / Frl * w
w3 <- fr / Frl * sqrt(w)
w4 <- fl / Frl * sqrt(w)
if (is.null(mmr))
mmr <- matrix(0, nrow = nrow(mml), ncol = ncol(mml))
if (is.null(mml))
mml <- matrix(0, nrow = nrow(mmr), ncol = ncol(mmr))
mmr[!is.finite(mmr)] <- 0
mml[!is.finite(mml)] <- 0
W3 <- mmr * w3
W4 <- mml * w4
return(-(crossprod(mmr * w1, mmr) - crossprod(mml * w2, mml) -
(crossprod(W3) - crossprod(W3, W4) -
crossprod(W4, W3) + crossprod(W4))))
}
.mlt_hessian_interval <- function(d, mml, mmr, offset = 0, mmtrunc = NULL,
w = 1) {
if (is.null(mmtrunc))
return(function(beta)
..mlt_hessian_interval(d, mml, mmr, offset, w, beta))
return(function(beta)
..mlt_hessian_interval(d, mml, mmr, offset, w, beta) -
..mlt_hessian_interval(d, mmtrunc$left,
mmtrunc$right, offset, w, beta))
}
.mlt_loglik_exact <- function(d, mm, mmprime, offset = 0, mmtrunc = NULL) {
if (is.null(mmtrunc))
return(function(beta)
d$d(offset + .xmb(mm, beta), log = TRUE) + .log(.xmb(mmprime, beta)))
return(function(beta)
d$d(offset + .xmb(mm, beta), log = TRUE) + .log(.xmb(mmprime, beta)) -
..mlt_loglik_interval(d, mmtrunc$left,
mmtrunc$right, offset, beta))
}
.mlt_score_exact <- function(d, mm, mmprime, offset = 0, mmtrunc = NULL) {
function(beta, Xmult = TRUE) {
mmb <- drop(.xmb(mm, beta)) + offset
ret <- d$dd2d(mmb)
if (Xmult)
ret <- ret * mm + mmprime / drop(.xmb(mmprime, beta))
if (!is.null(mmtrunc))
ret <- ret - ..mlt_score_interval(d, mmtrunc$left,
mmtrunc$right, offset, beta, Xmult)
return(ret)
}
}
.mlt_hessian_exact <- function(d, mm, mmprime, offset = 0, mmtrunc = NULL,
w = 1) {
function(beta) {
mmb <- drop(.xmb(mm, beta)) + offset
if (length(w) != length(mmb)) w <- rep(w, length(mmb))
w1 <- -(d$ddd(mmb) / d$d(mmb) - d$dd2d(mmb)^2) * w
w2 <- w / (drop(.xmb(mmprime, beta))^2)
ret <- crossprod(mm * w1, mm) + crossprod(mmprime * w2, mmprime)
if (!is.null(mmtrunc))
ret <- ret - ..mlt_hessian_interval(d, mmtrunc$left,
mmtrunc$right, offset, w, beta)
return(ret)
}
}
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.