Nothing
#' Fit generalized mixture cure model
#'
#' The following function fits a generalized mixture or non-mixture cure model
#' using a link function for the cure rate and for the survival of the uncured.
#' For a mixture cure model, the model is specified by
#' \deqn{S(t|z) = \pi(z) + [1 - \pi(z)] S_u(t|z),}
#' where
#' \deqn{g_1[S_u(t|z)] = \eta_1(t, z)\qquad and \qquad g_2[\pi(z)] = \eta_2(z).}
#' The function implements multiple link functions for both \eqn{g_1} and \eqn{g_2}. The default time-effect
#' is natural cubic splines, but the function allows for the use of other smoothers.
#'
#' @param formula Formula for modelling the survival of the uncured. Reponse has to be of the form \code{Surv(time, status)}.
#' @param data Data frame in which to interpret the variables names in \code{formula},
#' \code{smooth.formula}, \code{tvc.formula}, and \code{cr.formula}.
#' @param smooth.formula Formula for describing the time-effect of the survival of the uncured.
#' If \code{NULL} (default), the function uses the natural cubic splines modelled on the log-time scale.
#' @param smooth.args List. Optional arguments to the time-effect of the survival
#' of the uncured (default is \code{NULL}).
#' @param df Integer. Degrees of freedom (default is 3) for the time-effect of the survival of the uncured.
#' Is not used if \code{smooth.formula} is provided.
# @param logH.args
# @param logH.formula blabal
#' @param tvc Named list of integers. Specifies the degrees of freedom for time-varying covariate effects.
#' For instance, \code{tvc = list(x = 3)} creates a time-varying spline-effect of the covariate "x" with
#' 3 degrees of freedom using the \code{rstpm2::nsx} function.
#' @param tvc.formula Formula for the time-varying covariate effects.
#' For time-varying effects, a linear term of the covariate has to be included in \code{formula}.
#' @param cr.formula Formula for the cure proportion.
#' The left hand side of the formula is not used and should therefore not be specified.
#' @param bhazard Background hazard.
#' @param type A character indicating the type of cure model.
#' Possible values are \code{mixture} for mixture cure models (default) and \code{nmixture}
#' for non-mixture cure models.
#' @param covariance Logical. If \code{TRUE} (default), the covariance matrix is computed.
#' @param verbose Logical. If \code{TRUE} status messages of the function is outputted.
#' @param link.type.cr Character providing the link function for the cure proportion.
#' Possible values are \code{logit} (default), \code{loglog}, \code{identity}, and \code{probit}.
#' @param link.type Character providing the link function for the survival of the uncured.
#' Possible values are \code{PH} for a proportional hazards model (default), \code{PO} for a proportion odds model,
#' and \code{probit} for a probit model.
#' @param init Initial values for the optimization procedure.
#' If not specified, the function will create initial values internally.
#' @param timeVar Optional character giving the name of the variable specifying the time component of the \code{Surv} object.
#' Should currently not be used.
#' @param time0Var Optional character giving the name of the variable specifying the time start time component used for delayed entry.
#' Should currently not be used.
#' @param baseoff Logical. If \code{TRUE}, the time-effect is modelled only using \code{tvc.formula} rather
#' than merging with \code{smooth.formula}.
#' @param control Named list with control arguments passed to \code{optim}.
#' @param method Character passed to \code{optim} indicating the method for optimization.
#' See \code{?optim} for details.
#' @param constraint Logical. Indicates whether non-negativity constraints should be forced upon
#' the hazard of the uncured patients (see details).
#' @param ini.types Character vector denoting the executed schemes for computing initial values (see details).
# @param cure Logical. Indicates whether a cure model specification is needed for the survival of the uncured.
# This is usually \code{FALSE} (default).
#' @return An object of class \code{gfcm}.
#' @details The default smoother is natural cubic splines established by the \code{rstpm2::nsx} function.
#' Functions such as \code{ns}, \code{bs} are readily available for usage. Also the \code{cb} function in this package
#' can be used. Initial values are calculated by two procedures and the model is fitted under each set of initial values.
#' The model producing the highest likelihood is selected.\cr
#'
#' Using \code{link.type = 'PH'}, the link function \eqn{g_1(x) = \log(-\log(x))} is used.
#' Using \code{link.type = 'PO'}, the link function \eqn{g_1(x) = \log(\frac{x}{1 - x})} is used.\cr
#'
#' If \code{constraint = TRUE}, a non-negative hazard of the uncured is ensured by a general penalization scheme.
#' If \code{constraint = FALSE}, penalization is still employed, but on the all-cause hazard instead.
#' @export
#' @import survival
#' @import rstpm2
#' @import relsurv
#' @importFrom numDeriv hessian
#' @importFrom graphics abline
#' @importFrom stats .checkMFClasses as.formula delete.response dnorm lm model.frame model.matrix
#' na.pass optim pnorm printCoefmat pt qnorm quantile rnorm smooth.spline terms
#' @importFrom utils tail
#' @example inst/GenFlexCureModel.ex.R
#Note that the order of the terms in the smooth.formula matters.
#Meaning: smooth.formula = nsx(df = 4, x = something), will produce an error, while nsx(x = smothing, df= 4) will not
GenFlexCureModel <- function(formula, data, smooth.formula = NULL, smooth.args = NULL,
df = 3, tvc = NULL,
tvc.formula = NULL, bhazard = NULL, cr.formula = ~ 1,
type = "mixture",
link.type.cr = c("logit", "loglog", "identity", "probit"),
link.type = c("PH", "PO", "probit"),
init = NULL, baseoff = FALSE, timeVar = "", time0Var = "",
covariance = T, verbose = T,
control = list(maxit = 10000), method = "Nelder-Mead",
constraint = TRUE,
ini.types = c("cure", "flexpara")){
if(!type %in% c("mixture", "nmixture"))
stop("Wrong specication of argument 'type', must be either 'mixture' or 'nmixture'")
logH.args <- NULL
logH.formula <- NULL
cure <- FALSE
link.type <- match.arg(link.type)
link.surv <- switch(link.type, PH = link.PH, PO = link.PO, probit = link.probit,
AH = link.AH)
link.type.cr <- match.arg(link.type.cr)
if (!is.null(smooth.formula) && is.null(logH.formula))
logH.formula <- smooth.formula
if (!is.null(smooth.args) && is.null(logH.args))
logH.args <- smooth.args
eventInstance <- eval(lhs(formula),envir = data)
stopifnot(length(lhs(formula)) >= 2)
eventExpr <- lhs(formula)[[length(lhs(formula))]]
delayed <- length(lhs(formula)) >= 4
surv.type <- attr(eventInstance, "type")
if (surv.type %in% c("interval2", "left", "mstate"))
stop("stpm2 not implemented for Surv type ", surv.type,
".")
counting <- attr(eventInstance, "type") == "counting"
interval <- attr(eventInstance, "type") == "interval"
timeExpr <- lhs(formula)[[ifelse(delayed, 3, 2)]]
if (timeVar == "")
timeVar <- all.vars(timeExpr)
## set up the formulae
if (is.null(logH.formula) && is.null(logH.args)) {
logH.args$df <- df
if (cure)
logH.args$cure <- cure
}
if (is.null(logH.formula))
logH.formula <- as.formula(call("~", as.call(c(quote(nsx),
call("log", timeExpr), vector2call(logH.args)))))
if (is.null(tvc.formula) && !is.null(tvc)) {
tvc.formulas <- lapply(names(tvc), function(name) call(":",
as.name(name),
as.call(c(quote(nsx),
call("log",
timeExpr),
vector2call(if (cure) list(cure = cure,
df = tvc[[name]]) else list(df = tvc[[name]]))))))
if (length(tvc.formulas) > 1)
tvc.formulas <- list(Reduce(`%call+%`, tvc.formulas))
tvc.formula <- as.formula(call("~", tvc.formulas[[1]]))
}
if (!is.null(tvc.formula)) {
rhs(logH.formula) <- rhs(logH.formula) %call+% (rhs(tvc.formula))
}
if (baseoff)
rhs(logH.formula) <- rhs(tvc.formula)
full.formula <- formula
rhs(full.formula) <- rhs(formula) %call+% rhs(logH.formula)
.include <- apply(model.matrix(formula, data, na.action = na.pass),
1, function(row) !any(is.na(row))) & !is.na(eval(eventExpr,
data)) & !is.na(eval(timeExpr, data))
data <- data[.include, , drop = FALSE]
Call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "contrasts", "weights"),
names(mf), 0L)
mf <- mf[c(1L, m)]
time <- eval(timeExpr, data, parent.frame())
if (any(time > 0 & time < 1e-04))
warning("Some event times < 1e-4: consider transforming time to avoid problems with finite differences")
time0Expr <- NULL
if (delayed) {
time0Expr <- lhs(formula)[[2]]
if (time0Var == "")
time0Var <- all.vars(time0Expr)
time0 <- eval(time0Expr, data, parent.frame())
if (any(time0 > 0 & time0 < 1e-04))
warning("Some entry times < 1e-4: consider transforming time to avoid problems with finite differences")
}
event <- eval(eventExpr, data)
event <- if (length(unique(event)) == 1){
rep(TRUE, length(event))
} else {
event <- event > min(event)
}
# if (!interval) {
# coxph.call <- mf
# coxph.call[[1L]] <- as.name("coxph")
# coxph.call$model <- TRUE
# coxph.obj <- eval(coxph.call, envir = parent.frame())
# y <- model.extract(model.frame(coxph.obj), "response")
# data$logHhat <- pmax(-18, link.surv$link(Shat(coxph.obj)))
# }
# if (interval) {
# survreg.call <- mf
# survreg.call[[1L]] <- as.name("survreg")
# survreg.obj <- eval(survreg.call, envir = parent.frame())
# weibullShape <- 1/survreg.obj$scale
# weibullScale <- predict(survreg.obj)
# y <- model.extract(model.frame(survreg.obj), "response")
# data$logHhat <- pmax(-18, link$link(pweibull(time, weibullShape,
# weibullScale, lower.tail = FALSE)))
# }
# lm.call <- mf
# lm.call[[1L]] <- as.name("lm")
# lm.formula <- full.formula
# lhs(lm.formula) <- quote(logHhat)
# lm.call$formula <- lm.formula
# dataEvents <- data[event, ]
# if (interval)
# dataEvents <- data
# lm.call$data <- quote(dataEvents)
# lm.obj <- eval(lm.call)
# mt <- terms(lm.obj)
# mf <- model.frame(lm.obj)
lm.call <- mf
lm.call[[1L]] <- as.name("lm")
lm.formula <- full.formula
lhs(lm.formula) <- quote(arbri)
lm.call$formula <- lm.formula
dataEvents <- data[event, ]
dataEvents$arbri <- rnorm(nrow(dataEvents))
if (interval)
dataEvents <- data
lm.call$data <- quote(dataEvents)
lm.obj <- eval(lm.call)
#Create background hazard
if(is.null(bhazard)){
bhazard <- rep(0, nrow(data))
}else {
if(!is.numeric(bhazard)){
bhazard <- data[, bhazard]
}
}
excess <- !all(bhazard == 0)
if(length(bhazard) != nrow(data))
stop("Length of bhazard is not the same as nrow(data)")
lpfunc <- function(delta, fit, dataset, var) {
dataset[[var]] <- dataset[[var]] + delta
lpmatrix.lm(fit, dataset)
}
## initialise values specific to either delayed entry or interval-censored
ind0 <- FALSE
map0 <- 0L
which0 <- 0
#wt0 <- 0
ttype <- 0
transX <- function(X, data) X
transXD <- function(XD) XD
if (!interval) {
X <- lpmatrix.lm(lm.obj, data)
if (link.type == "AH") {
datat0 <- data
datat0[[timeVar]] <- 0
index0 <- which.dim(X - lpmatrix.lm(lm.obj, datat0))
transX <- function(X, data) {
datat0 <- data
datat0[[timeVar]] <- 0
Xt0 <- lpmatrix.lm(lm.obj, datat0)
(X - Xt0)[, index0, drop = FALSE]
}
transXD <- function(XD) XD[, index0, drop = FALSE]
init <- init[index0]
}
X <- transX(X, data)
XD <- grad(lpfunc, 0, lm.obj, data, timeVar)
XD <- transXD(matrix(XD, nrow = nrow(X)))
X1 <- matrix(0, nrow(X), ncol(X))
X0 <- matrix(0, 1, ncol(X))
if (delayed && all(time0 == 0))
delayed <- FALSE
if (delayed) {
ind0 <- time0 > 0
map0 <- vector("integer", nrow(X))
map0[ind0] <- as.integer(1:sum(ind0))
map0[!ind0] <- NaN
which0 <- 1:nrow(X)
which0[!ind0] <- NaN
data0 <- data[ind0, , drop = FALSE]
data0[[timeVar]] <- data0[[time0Var]]
X0 <- transX(lpmatrix.lm(lm.obj, data0), data0)
#wt0 <- wt[ind0]
rm(data0)
}
} else {
ttype <- eventInstance[, 3]
X1 <- transX(lpmatrix.lm(lm.obj, data), data)
data0 <- data
data0[[timeVar]] <- data0[[time0Var]]
X <- transX(lpmatrix.lm(lm.obj, data0), data0)
XD <- grad(lpfunc, 0, lm.obj, data0, timeVar)
XD <- transXD(matrix(XD, nrow = nrow(X)))
X0 <- matrix(0, nrow(X), ncol(X))
rm(data0)
}
#Create linear object for the cure rate part
lm.call <- mf
lm.call[[1L]] <- as.name("lm")
lm.formula <- cr.formula
lhs(lm.formula) <- quote(arbri)
lm.call$formula <- lm.formula
dataEvents <- data
dataEvents$arbri <- rnorm(nrow(dataEvents))
if (interval)
dataEvents <- data
lm.call$data <- quote(dataEvents)
lm.obj.cr <- eval(lm.call)
X.cr <- lpmatrix.lm(lm.obj.cr, data)
#X.cr <- model.matrix(cr.formula, data = data)
if(is.null(init)){
if(verbose) cat("Finding initial values... ")
#ini.types <- if(delayed) ini.types[2] else ini.types
init <- vector("list", length(ini.types))
for(i in 1:length(init)){
args <- list(formula = formula, data = data, smooth.formula = smooth.formula,
logH.formula = logH.formula, tvc.formula = tvc.formula, cr.formula = cr.formula,
full.formula = full.formula, X = X, X0 = X0, X.cr = X.cr, delayed = delayed,
bhazard = bhazard, type = type, link.type.cr = link.type.cr,
link.surv = link.surv, time = time, timeExpr = as.character(timeExpr),
lm.obj = lm.obj, method = ini.types[i])
init[[i]] <- do.call(get.init, args)
}
} else {
if(verbose) cat("Initial values provided by the user... ")
}
# #Extract minus log likelihood function
# if(delayed) {
# minusloglik <- switch(type,
# mixture = GenFlexMixMinLogLikDelayed,
# nmixture = GenFlexNmixMinLogLikDelayed)
# } else {
# minusloglik <- switch(type,
# mixture = GenFlexMixMinLogLik,
# nmixture = GenFlexNmixMinLogLik)
# }
cure.type <- switch(type,
mixture = mix,
nmixture = nmix)
#minusloglik <- switch(type,
# mixture = GenFlexMixMinLogLikDelayed,
# nmixture = GenFlexNmixMinLogLikDelayed)
minusloglik <- GenFlexMinLogLikDelayed
#Prepare optimization arguments
args <- list(event = event, X = X, XD = XD, X.cr = X.cr, X0 = X0, ind0 = ind0,
bhazard = bhazard, link.type.cr = link.type.cr,
link.surv = link.surv, kappa = 0, constraint = FALSE, cure.type = cure.type)
if(is.null(control$maxit)){
control$maxit <- 10000
}
#Test if initial values are within the feasible region
ini.eval <- sapply(init, function(inival) do.call(minusloglik, c(args, list(inival))))
run.these <- !is.na(ini.eval)
if(all(!run.these))
stop("Initial values are outside feasible region")
if(verbose) cat("Completed!\nFitting the model... ")
optim.args <- c(control = list(control), args)
optim.args$kappa <- 1
optim.args$fn <- minusloglik
optim.args$method <- method
optim.args$constraint <- constraint
optim.args$cure.type <- cure.type
res_list <- vector("list", length(init[run.these]))
for(i in 1:length(res_list)){
neghaz <- T
while(neghaz){
#cat(optim.args$kappa, "\n")
optim.args$par <- init[run.these][[i]]
res.optim <- do.call(optim, optim.args)
beta <- res.optim$par[(ncol(X.cr) + 1):length(res.optim$par)]
eta <- X %*% beta
etaD <- XD %*% beta
if(constraint){
haz.const <- link.surv$h(eta, etaD)
} else {
gamma <- res.optim$par[1:ncol(X.cr)]
eta.pi <- X.cr %*% gamma
pi <- get.link(link.type.cr)(eta.pi)
surv <- link.surv$ilink(eta)
rsurv <- cure.type$surv(pi, surv)
ehaz <- cure.type$haz(pi, link.surv$gradS(eta, etaD), rsurv)
haz.const <- bhazard + ehaz
}
neghaz <- any(haz.const < 0)
optim.args$kappa <- optim.args$kappa * 10
}
optim.args$kappa <- 1
res_list[[i]] <- res.optim
}
#Choose the best model according to the maximum likelihood estimate
MLs <- sapply(res_list, function(x) tail(x$value, 1))
wh <- which.min(MLs)
res <- res_list[[wh]]
if(res$convergence != 0){
warning("Convergence not reached")
}
if(verbose) cat("Completed!\n")
#Compute the covariance matrix matrix
if(covariance){
args$kappa <- 0
args$x <- res$par
args$func <- minusloglik
args$constraint <- FALSE
hes <- do.call(numDeriv::hessian, args)
cov <- if (!inherits(vcov <- try(solve(hes)), "try-error")) vcov
# cov <- solve(hes)
if(!is.null(cov) && any(is.na(cov))){
warning("Hessian is not invertible!")
}
}else{
cov <- NULL
}
#Output the results
L <- list(formula = formula, smooth.formula = smooth.formula, tvc.formula = tvc.formula,
logH.formula = logH.formula, cr.formula = cr.formula, full.formula = full.formula,
coefs = res$par[1:ncol(X.cr)],
coefs.spline = res$par[(ncol(X.cr) + 1):length(res$par)],
data = data, NegMaxLik = min(MLs), covariance = cov, ci = covariance,
type = type, NegMaxLiks = MLs, optim.pars = optim.args[c("control", "fn")],
args = args, timeExpr = timeExpr, lm.obj = lm.obj, lm.obj.cr = lm.obj.cr,
link.type.cr = link.type.cr, link.type = link.type, link.surv = link.surv, excess = excess,
timeVar = timeVar, transX = transX, transXD = transXD,
time = time, event = event, eventExpr = eventExpr, cure.type = cure.type, ML = -min(MLs))
class(L) <- c("gfcm", "cuRe")
L
}
#Function for computing initial values
get.init <- function(formula, data, smooth.formula, logH.formula, tvc.formula, cr.formula, full.formula, delayed,
bhazard, type, link.type.cr, link.surv, timeExpr, time, lm.obj, X, X0, X.cr, method){
if(!method %in% c("cure", "flexpara")){
stop("Argument method should be either 'cure' or 'flexpara'")
}
if(method == "cure"){
formula.pi <- cr.formula
lhs(formula.pi) <- lhs(formula)
formula.k1 <- formula
lhs(formula.k1) <- NULL
if(length(attr(terms(formula.k1), "term.labels"))){
a <- Reduce(paste, deparse(formula.k1))
a <- gsub("-1", "1", a)
formula.k1 <- as.formula(a)
} else {
formula.k1 <- ~ 1
}
#Fit mixture or non-mixture cure model
fit <- fit.cure.model(formula = formula.pi, data = data, bhazard = bhazard, covariance = F,
formula.surv = list(formula.k1, ~ 1), type = type)
#Scale by link function
pi_hat <- do.call(rbind, predict(fit, type = "curerate", newdata = data))$Estimate
#Predict survival of the uncured
lp <- exp(model.matrix(formula.k1, data = data) %*% fit$coefs[[2]])
suhat <- exp(-lp * time ^ exp(fit$coefs[[3]]))
}else if(method == "flexpara"){
formula.2 <- formula
vars1 <- attr(terms(cr.formula), "term.labels")
vars2 <- attr(terms(formula.2), "term.labels")
wh <- which(!vars1 %in% vars2)
if(length(wh)){
formula.pi <- as.formula(paste0("~ ", paste(vars1, collapse = "+")))
rhs(formula.2) <- rhs(formula) %call+% rhs(formula.pi)
}
if(length(attr(terms(formula.2), "term.labels"))){
a <- Reduce(paste, deparse(formula.2))
a <- gsub("-1", "1", a)
formula.2 <- as.formula(a)
} else {
rhs(formula.2) <- 1
}
#Fit relative survival model
suppressWarnings(fit <- do.call(rstpm2::stpm2, list(formula = formula.2, data = data, bhazard = bhazard)))
#Predict survival function
shat <- predict(fit, newdata = data, se.fit = F, keep.attributes = F)
#If predictions are all 1, we manually change these
shat[shat == 1] <- shat[shat == 1] - 0.01
#Change follow-up times and predict cure rate
data2 <- data
data2[, timeExpr] <- max(data2[, timeExpr]) + 0.1
pi_hat <- predict(fit, newdata = data2, se.fit = F, keep.attributes = F)
#Change cases with increasing relative survival
wh <- which(pi_hat >= shat)
pi_hat[wh] <- shat[wh] - 0.01
#Run linear model for S_u(t) to obtain initial values for either mixture or non-mixture models
if(type == "mixture"){
suhat <- (shat - pi_hat) / (1 - pi_hat)
} else {
suhat <- 1 - log(shat) / log(pi_hat)
}
}
#Run linear model for pi to obtain initial values
gpi_hat <- get.inv.link(link.type.cr)(pi_hat)
ini_pi <- lm(gpi_hat ~ -1 + X.cr)$coefficients
names(ini_pi) <- colnames(X.cr)
#Run linear model for survival to obtain initial values
gsuhat <- link.surv$link(suhat)
finites <- is.finite(gsuhat)
suppressWarnings(ini_surv <- lm(gsuhat[finites] ~ -1 + X[finites,])$coefficients)
names(ini_surv) <- colnames(X)
c(ini_pi, ini_surv)
}
#' @export
#' @method print gfcm
#Print function for class gfcm
print.gfcm <- function(x, ...){
cat("Call pi:\n")
print(x$formula)
cat("Call S_u(t):\n")
print(x$formula_main)
cat("\nCoefficients:\n")
print(list(pi = x$coefs,
surv = x$coefs.spline))
}
#' @export
#' @method summary gfcm
#Summary function for class gfcm
summary.gfcm <- function(object, ...){
se <- sqrt(diag(object$covariance))
zval <- c(object$coefs, object$coefs.spline) / se
TAB <- cbind(Estimate = c(object$coefs, object$coefs.spline),
StdErr = se,
z.value = zval,
p.value = ifelse(is.na(zval), rep(NA, length(se)),
2 * pnorm(-abs(zval))))
TAB1 <- TAB[1:length(object$coefs), ,drop = F]
TAB2 <- TAB[(length(object$coefs) + 1):(length(object$coefs.spline) + length(object$coefs)),, drop = F]
results <- list(pi = TAB1, surv = TAB2)
results$type <- object$type
results$linkpi <- object$link.type.cr
results$linksu <- object$link.type
results$ML <- object$NegMaxLik
results$formula <- object$formula
results$smooth.formula <- object$smooth.formula
results$cr.formula <- object$cr.formula
results$tvc.formula <- object$tvc.formula
results$full.formula <- object$full.formula
class(results) <- "summary.gfcm"
results
}
#' @export
#' @method print summary.gfcm
#Print for class summary.gfcm
print.summary.gfcm <- function(x, ...)
{
cat("Call - pi:\n")
print(x$formula)
# cat("\n")
stats::printCoefmat(x$pi, P.values = TRUE, has.Pvalue = T, signif.legend = F)
cat("\nCall - surv:\n")
print(x$full.formula)
# if(length(all.vars(x$formula.tvc))){
# cat("Call - surv - tvc: ")
# print(deparse(x$formula.tvc))
# }
stats::printCoefmat(x$surv, P.values = TRUE, has.Pvalue = T)
cat("\n")
cat("Type =", x$type, "\n")
cat("Link - pi =", x$linkpi, "\n")
cat("Link - surv = ", x$linksu, "\n")
cat("LogLik(model) =", x$ML, "\n")
}
#Functionalities from the rstpm2 package - thanks to Mark Clements for brilliant code
## link families
link.PH <- list(link=function(S) log(-log(as.vector(S))),
ilink=function(eta) exp(-exp(as.vector(eta))),
gradS=function(eta,X) -exp(as.vector(eta))*exp(-exp(as.vector(eta)))*X,
h=function(eta,etaD) as.vector(etaD)*exp(as.vector(eta)),
H=function(eta) exp(as.vector(eta)),
gradh=function(eta,etaD,obj) obj$XD*exp(as.vector(eta))+obj$X*as.vector(etaD)*exp(as.vector(eta)),
gradH=function(eta,obj) obj$X*exp(as.vector(eta)))
link.PO <- list(link=function(S) -logit(as.vector(S)),
ilink=function(eta) expit(-as.vector(eta)),
gradS=function(eta,X) -(exp(as.vector(eta))/(1+exp(as.vector(eta)))^2)*X,
H=function(eta) log(1+exp(as.vector(eta))),
h=function(eta,etaD) as.vector(etaD)*exp(as.vector(eta))*expit(-as.vector(eta)),
gradh=function(eta,etaD,obj) {
as.vector(etaD)*exp(as.vector(eta))*obj$X*expit(-as.vector(eta)) -
exp(2*as.vector(eta))*obj$X*as.vector(etaD)*expit(-as.vector(eta))^2 +
exp(as.vector(eta))*obj$XD*expit(-as.vector(eta))
},
gradH=function(eta,obj) obj$X*exp(as.vector(eta))*expit(-as.vector(eta)))
link.probit <-
list(link=function(S) -qnorm(as.vector(S)),
ilink=function(eta) pnorm(-as.vector(eta)),
gradS=function(eta,X) -dnorm(-as.vector(eta))*X,
H=function(eta) -log(pnorm(-as.vector(eta))),
h=function(eta,etaD) dnorm(as.vector(eta))/pnorm(-as.vector(eta))*as.vector(etaD),
gradh=function(eta,etaD,obj) {
-as.vector(eta)*obj$X*dnorm(as.vector(eta))*as.vector(etaD)/pnorm(-as.vector(eta)) +
obj$X*dnorm(as.vector(eta))^2/pnorm(-as.vector(eta))^2*as.vector(etaD) +
dnorm(as.vector(eta))/pnorm(-as.vector(eta))*obj$XD
},
gradH=function(eta,obj) obj$X*dnorm(as.vector(eta))/pnorm(-as.vector(eta)))
link.AH <- list(link=function(S) -log(S),
ilink=function(eta) exp(-as.vector(eta)),
gradS=function(eta,X) -as.vector(exp(-as.vector(eta)))*X,
h=function(eta,etaD) as.vector(etaD),
H=function(eta) as.vector(eta),
gradh=function(eta,etaD,obj) obj$XD,
gradH=function(eta,obj) obj$X)
rhs=function(formula)
if (length(formula)==3) formula[[3]] else formula[[2]]
lhs <- function(formula)
if (length(formula)==3) formula[[2]] else NULL
"rhs<-" = function(formula,value) {
newformula <- formula
newformula[[length(formula)]] <- value
newformula
}
"lhs<-" <- function(formula,value) {
if (length(formula)==2)
as.formula(as.call(c(formula[[1]],value,formula[[2]])))
else {
newformula <- formula
newformula[[2]] <- value
newformula
}
}
allCall=function(obj) {
if (is.atomic(obj) && length(obj)==1) return(obj)
if (is.atomic(obj) && length(obj)>1) return(as.call(c(quote(c),as.list(obj))))
if (is.name(obj) || is.symbol(obj)) return(obj)
as.call(lapply(obj,allCall))
}
vector2call=function(obj) {
if (is.atomic(obj) && length(obj)==1) return(obj)
if (is.atomic(obj) && length(obj)>1) return(as.call(c(quote(c),as.list(obj))))
if (is.name(obj) || is.symbol(obj)) return(obj)
lapply(obj,allCall) # is this correct?
}
`%call+%` <- function(left,right) call("+",left,right)
## predict lpmatrix for an lm object
lpmatrix.lm <-
function (object, newdata, na.action = na.pass) {
tt <- terms(object)
if (!inherits(object, "lm"))
warning("calling predict.lm(<fake-lm-object>) ...")
if (missing(newdata) || is.null(newdata)) {
X <- model.matrix(object)
}
else {
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.action,
xlev = object$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses")))
.checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
}
X
}
## numerically calculate the partial gradient \partial func_j \over \partial x_i
## (dim(grad(func,x)) == c(length(x),length(func(x)))
grad <- function(func,x,...) # would shadow numDeriv::grad()
{
h <- .Machine$double.eps^(1/3)*ifelse(abs(x)>1,abs(x),1)
temp <- x+h
h.hi <- temp-x
temp <- x-h
h.lo <- x-temp
twoeps <- h.hi+h.lo
nx <- length(x)
ny <- length(func(x,...))
if (ny==0L) stop("Length of function equals 0")
df <- if(ny==1L) rep(NA, nx) else matrix(NA, nrow=nx,ncol=ny)
for (i in 1L:nx) {
hi <- lo <- x
hi[i] <- x[i] + h.hi[i]
lo[i] <- x[i] - h.lo[i]
if (ny==1L)
df[i] <- (func(hi, ...) - func(lo, ...))/twoeps[i]
else df[i,] <- (func(hi, ...) - func(lo, ...))/twoeps[i]
}
return(df)
}
which.dim <- function (X, silent = TRUE)
{
stopifnot(is.matrix(X))
silent <- as.logical(silent)[1]
qr.X <- qr(X, tol = 1e-07, LAPACK = FALSE)
if (qr.X$rank == ncol(X))
return(TRUE)
if (!silent)
message(gettextf("design is column rank deficient so dropping %d coef",
ncol(X) - qr.X$rank))
return(qr.X$pivot[1:qr.X$rank])
}
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.