R/CoxRFX.R
In gerstung-lab/ebmstate: Empirical Bayes Multi-State Cox Model
Defines functions
strata
MakeInteger
print.msfit
print.coxrfx
summary.coxrfx
CoxRFX
Documented in
CoxRFX
MakeInteger
print.coxrfx
print.msfit
summary.coxrfx
#' Empirical Bayes, multi-state Cox model
#'
#' This function estimates a multi-state Cox model with one or more Gaussian priors
#' imposed on the regression coefficients (see Therneau et al., 2003).
#' Multiple groups of coefficients can be defined: coefficients within a group share
#' the same (possibly unknown) mean and variance. The parameters and hyperparameters are
#' efficiently estimated by an EM-type algorithm built around the function
#' \code{survival::coxph}.
#' @param Z A data frame consisting of the covariate columns of a data set in 'long format',
#' and two extra columns: one named 'trans', with
#' the transition that each row refers to, and
#' another named 'strata', with the stratum
#' of each transition (transitions belonging to the same
#' stratum are assumed to have the same baseline hazard
#' function).
#' @param surv A `survival' object created with \code{survival::Surv}.
#' @param groups A character or numeric vector whose \eqn{i}th element gives the group of the regression
#' coefficient associated with the \eqn{i}th covariate column of Z (coefficients belonging to the same group
#' share the same Gaussian prior).
#' @param tmat Transition matrix describing the states and transitions in
#' the multi-state model. See \code{trans} in \code{\link[mstate:msprep]{mstate::msprep}}
#' for more detailed information.
#' @param which.mu A vector with names or numbers of coefficient groups (see
#' argument \code{groups}). If the name or number of a group of coefficients is
#' given in this argument, \code{CoxRFX} will estimate the mean of its Gaussian distribution;
#' otherwise the mean will be fixed at zero.
#' @param tol Convergence criterium of the EM algorithm. The algorithm stops unless
#' there is at least one parameter (or hyperparameter) for which it holds that the
#' current estimate differs in absolute terms by more than \code{tol} from the
#' previous estimate.
#' @param max.iter The maximum number of iterations in the EM algorithm.
#' @param sigma0 A vector with the initial value of the variance hyperparameter for each group of coefficients.
#' Or a single value, in case the initial value of the variance hyperparameter is meant to be the same for all groups.
#' @param sigma.hat Which estimator to use for the variance hyperparameters (see details).
#' @param verbose Gives more output.
#' @param ... Further arguments passed to the function \code{survival::coxph}.
#' @details
#'
#' Different estimators exist for the variance hyperparameters: the default is "df", as used by Perperoglou (2014) and introduced by Schall (1991).
#' Alternatives are MLE, REML, and BLUP, as defined by Therneau et al. (2003).
#' Simulations suggest that the 'df' method is the most accurate.
#'
#' The model can also be fitted using package \code{coxme}; the \code{coxme}
#' routine numerically optimises the integrated partial likelihood, which may
#' be more accurate, but is computationally expensive.
#'
#' @references Terry M Therneau, Patricia M Grambsch & V. Shane Pankratz (2003) Penalized Survival Models and Frailty, Journal of Computational and Graphical Statistics, 12:1, 156-175, http://dx.doi.org/10.1198/1061860031365
#'
#' A. Perperoglou (2014). Cox models with dynamic ridge penalties on time-varying effects of the covariates. Stat Med, 33:170-80. http://dx.doi.org/10.1002/sim.5921
#'
#' R. Schall (1991). Estimation in generalized linear models with random effects. Biometrika, 78:719-727. http://dx.doi.org/10.1093/biomet/78.4.719
#' @return An object of class \code{c(coxrfx,coxph.penal,coxph)},
#' which is essentially a \code{coxph} object with a few extra
#' fields [the inputs $groups, $Z, $surv, and $tmat,
#' and the hyperparameters $sigma2 (variances) and $mu (means)].
#' See \code{\link[survival:coxph.object]{survival::coxph.object}}.
#'
#' @author Moritz Gerstung & Rui Costa, extending the work of
#' Terry Therneau et al. in the package \code{survival}.
#' @seealso Package \code{survival}
#' \code{\link[survival:coxph.object]{survival::coxph.object}};
#' \code{\link[survival:Surv]{survival::Surv}}; package \code{coxme}.
#' @export
#' @example inst/examples/CoxRFX-example.R
CoxRFX <- function(Z, surv, groups = rep(1, ncol(Z)),tmat=NULL, which.mu = unique(groups), tol=1e-3, max.iter=50, sigma0 = 0.1, sigma.hat=c("df","MLE","REML","BLUP"), verbose=FALSE, ...){
##
trans<-Z$trans
strata<-Z$strata
Z$strata<-Z$trans<-NULL
namesZ<-names(Z)
##
Z = as.matrix(Z)
Z.df <- TRUE
if(is.null(colnames(Z)))
colnames(Z) <- make.names(1:ncol(Z))
sigma.hat = match.arg(sigma.hat)
o <- order(groups)
Z <- Z[,o]
groups <- factor(groups[o])
uniqueGroups <- levels(groups)
ZZ <- lapply(uniqueGroups, function(i) Z[,groups==i, drop=FALSE])
names(ZZ) <- uniqueGroups
sumZ <- sapply(which.mu, function(i) rowSums(ZZ[[i]]))
nGroups = length(uniqueGroups)
sigma2 <- sigma0ld <- rep(ifelse(sigma0>0, sigma0,1), nGroups)
iter = 1
mu <- mu0ld <- rep(0, nGroups)
names(mu) <- uniqueGroups
beta = rep(1,ncol(Z)+length(which.mu))
beta0ld = rep(0,ncol(Z)+length(which.mu))
sigma2.mu = sigma0
nu=0
penalize.mu = FALSE
if(!is.null(which.mu))
if(!penalize.mu)
sumTerm <- "sumZ"
else
sumTerm <- "ridge(sumZ, theta=1/sigma2.mu, scale=FALSE)"
else sumTerm <- character(0)
while((max(abs(beta-beta0ld)) > tol | max(abs(mu - mu0ld)) > tol | max(abs(sigma2 - sigma0ld)) > tol) & iter < max.iter){
beta0ld = beta
sigma0ld <- sigma2
mu0ld <- mu
formula <- formula(paste("surv ~", paste(c(sapply(1:nGroups, function(i) paste("ridge(ZZ[[",i,"]], theta=1/sigma2[",i,"], scale=FALSE)", sep="")),
sumTerm,"strata(strata)"),
collapse=" + ")))
fit <- coxph(formula, ...)
if(any(is.na(coef(fit)))){
warning(paste("NA during estimation (iter: ", iter, ", coef: ", paste(which(is.na(coef(fit)[order(o)])), sep=","), ")", sep=""))
break
}
if(!is.null(which.mu))
mu[which.mu] <- coef(fit)[-(1:ncol(Z))]
if(verbose) cat("mu", mu, "\n", sep="\t")
names(fit$df) <- c(uniqueGroups, rep("Offset", length(which.mu)>0))
if(verbose) cat("df", fit$df,"\n", sep="\t")
sigma2 = sapply(uniqueGroups, function(i){
index <- which(groups==i) #& fit$coefficients > beta.thresh
if(sigma.hat=="BLUP")
(nu * sigma0 + sum((fit$coefficients[index])^2 ))/(nu + length(index))
else if(sigma.hat=="df")
(nu * sigma0 + sum((fit$coefficients[index])^2 ))/(nu + fit$df[i])
else if(sigma.hat == "MLE")
(nu * sigma0 + sum((fit$coefficients[index])^2 ) + sum(diag(solve(solve(fit$var)[index,index]))))/(nu + length(index))
else if(sigma.hat == "REML")
(nu * sigma0 + sum((fit$coefficients[index])^2 ) + sum(diag(fit$var)[index]))/(nu + length(index))
})
if(verbose) {
cat("sigma2", sigma2, "\n", sep="\t")
cat("loglik:", fit$loglik - c(0,fit$penalty[2] + 1/2 * sum(log(sigma2[groups]))),"\n", sep="\t")
}
if(penalize.mu){
if(sigma.hat=="BLUP")
sigma2.mu = (sigma0 * nu + sum((mu-0)^2)) / (nu + length(mu))
else if(sigma.hat=="df")
sigma2.mu = (sigma0 * nu + sum((mu-0)^2)) / (nu + fit$df["Offset"])
else if(sigma.hat == "MLE")
sigma2.mu = (nu * sigma0 + sum((mu-0)^2 ) + sum(diag(solve(solve(fit$var)[-(1:ncol(Z)),-(1:ncol(Z))]))))/(nu + length(mu))
else if(sigma.hat == "REML")
sigma2.mu = (nu * sigma0 + sum((mu-0)^2 ) + sum(diag(fit$var)[-(1:ncol(Z))]))/(nu + length(mu))
}
beta = fit$coefficients
iter = iter+1
}
if(iter == max.iter)
warning("Did not converge after ", max.iter, " iterations.")
fit$iter[1] <- iter
fit$sigma2 = sigma0ld
names(fit$sigma2) <- uniqueGroups
#fit$sigma2.mu = sigma2.mu
fit$mu = mu
fit$Z = data.frame(Z[,order(o)],trans=trans,strata=strata)
fit$surv = surv
fit$tmat=tmat
C <- rbind(diag(1, ncol(Z)),t(as.matrix(MakeInteger(groups)[which.mu]))) ## map from centred to uncentred coefficients
fit$groups = groups[order(o)]
var = fit$var
var2 = fit$var2
colnames(var) <- rownames(var) <- colnames(var2) <- rownames(var2) <- rownames(C) <- c(colnames(Z), which.mu)
colnames(C) <- colnames(Z)
p <- ncol(Z)
i <- c(order(o), (1:ncol(var))[-p:-1])
j <- order(o)
fit$C <- C[i,j]
fit$Hinv <- var[i,i] ## Hinv
fit$V <- var2[i,i] ## Hinv I Hinv
fit$z <- (fit$coefficients / sqrt(diag(var)))[i] ## z-scores of centred coefficients
fit$z2 <- (fit$coefficients / sqrt(diag(var2)))[i] ## z-scores of centred coefficients (var2)
fit$var = (t(C) %*% var %*% C)[j,j] ## covariance of uncentred coef
fit$var2 = (t(C) %*% var2 %*% C)[j,j] ## covariance of uncentred coef (var2)
fit$mu.var = var[-(1:p),-(1:p)] ## covariance of mean
fit$mu.var2 = var2[-(1:p),-(1:p)] ## covariance of mean (var2)
fit$means = fit$means[1:p][j]
fit$coefficients <- (fit$coefficients %*% C)[j]
names(fit$means) <- names(fit$coefficients) <- namesZ
fit$terms <- fit$terms[1:length(uniqueGroups)]
fit$penalized.loglik <- fit$loglik[2] - fit$penalty[2] - 1/2 * sum(log(fit$sigma2[groups]))
## Fake call for predict.coxph and survfit.coxph
call <- match.call()
if(Z.df){
call["data"] <- call["Z"]
formula <- as.formula(paste(as.character(call["surv"]),"~",paste(colnames(Z)[j], collapse="+"),"+strata(strata)"))
}else{
formula <- as.formula(paste(as.character(call["surv"]),"~",as.character(call["Z"])))
}
attr(formula,".Environment") <- parent.frame()
fit$formula <- formula
call["formula"] <- call("foo",formula=formula)["formula"]
fit$terms <- terms(formula)
attr(fit$terms,"specials")<-list(strata=NULL,cluster=NULL)
attr(fit$terms,"specials")$strata<-length(attr(fit$terms,"term.labels"))+1
fit$call <- call
class(fit) <- c("coxrfx", class(fit))
return(fit)
}
#' A summary method for CoxRFX models
#'
#' This function prints the point estimates of parameters and
#' hyperparameters contained in a \code{coxrfx} object.
#'
#'
#' @param object A \code{coxrfx} object
#' (obtained by running the function \code{CoxRFX}).
#' @param ... Further arguments passed to or from other methods.
#' @return NULL
#'
#' @author Rui Costa
#' @export
#' @method summary coxrfx
summary.coxrfx <- function(object,...){
cat("Hyperparameters:\n")
show(format(data.frame(mean=object$mu, variance=object$sigma2), digits=2))
cat("\n Parameters: \n")
show(format(data.frame(group=object$groups, estimate=object$coefficients,row.names = names(object$coefficients)), digits=2))
}
#' Print method for CoxRFX objects
#'
#' This function implicitly calls summary.coxrfx().
#' @param x A \code{coxrfx} object
#' @param ... further arguments passed to or from other methods.
#' @return NULL
#'
#' @author Moritz Gerstung & Rui Costa
#' @export
#' @method print coxrfx
print.coxrfx <- function(x,...){
summary.coxrfx(x)
}
#' Print method for \code{msfit} objects
#' generated by \code{msfit_generic}
#'
#' This method is a simple call to \code{print.default}.
#' Its main purpose is to override \code{print.coxrfx}
#' when printing an
#' object of double class \code{msfit} and \code{coxrfx}.
#' @param x An object of class \code{msfit} or double class \code{msfit}
#' and \code{coxrfx}.
#' @param ... Further arguments passed to or from other methods.
#' @return NULL
#'
#' @author Rui Costa
#' @export
#' @method print msfit
print.msfit <- function(x,...){
#mstate:::summary.msfit(x)
print.default(x,...)
}
#' Convert factor to integer.
#' @param v A factor vector.
#' @return A data.frame with columns corresponding to levels in the factor.
#' @details An internal function of \code{CoxRFX}, not meant
#' to called directly by the user.
#' @author Moritz Gerstung
#' @seealso \code{\link{CoxRFX}}
MakeInteger <- function(v){
res <- as.data.frame(lapply(levels(v), `==`, v))
colnames(res) <- levels(v)
res + 0
}
#coxph function from survival package version 2.44-1.1
#followed by non-exported functions of the same package (potentially used by coxph)
coxph<-function (formula, data, weights, subset, na.action, init, control,
ties = c("efron", "breslow", "exact"), singular.ok = TRUE,
robust = FALSE, model = FALSE, x = FALSE, y = TRUE, tt, method = ties,
...)
{
ties <- match.arg(ties)
Call <- match.call()
extraArgs <- list(...)
if (length(extraArgs)) {
controlargs <- names(formals(coxph.control))
indx <- pmatch(names(extraArgs), controlargs, nomatch = 0L)
if (any(indx == 0L))
stop(gettextf("Argument %s not matched", names(extraArgs)[indx ==
0L]), domain = NA)
}
if (missing(control))
control <- coxph.control(...)
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
names(Call), nomatch = 0)
if (indx[1] == 0)
stop("A formula argument is required")
temp <- Call[c(1, indx)]
temp[[1L]] <- quote(stats::model.frame)
special <- c("strata", "cluster", "tt")
temp$formula <- if (missing(data))
terms(formula, special)
else terms(formula, special, data = data)
if (!is.null(attr(temp$formula, "specials")$tt)) {
coxenv <- new.env(parent = environment(formula))
assign("tt", function(x) x, env = coxenv)
environment(temp$formula) <- coxenv
}
mf <- eval(temp, parent.frame())
if (nrow(mf) == 0)
stop("No (non-missing) observations")
Terms <- terms(mf)
Y <- model.extract(mf, "response")
if (!inherits(Y, "Surv"))
stop("Response must be a survival object")
type <- attr(Y, "type")
if (type != "right" && type != "counting")
stop(paste("Cox model doesn't support \"", type, "\" survival data",
sep = ""))
data.n <- nrow(Y)
if (control$timefix)
Y <- aeqSurv(Y)
if (length(attr(Terms, "variables")) > 2) {
ytemp <- terms.inner(formula[1:2])
xtemp <- terms.inner(formula[-2])
if (any(!is.na(match(xtemp, ytemp))))
warning("a variable appears on both the left and right sides of the formula")
}
strats <- attr(Terms, "specials")$strata
if (length(strats)) {
stemp <- untangle.specials(Terms, "strata", 1)
if (length(stemp$vars) == 1)
strata.keep <- mf[[stemp$vars]]
else strata.keep <- strata(mf[, stemp$vars], shortlabel = TRUE)
strats <- as.numeric(strata.keep)
}
timetrans <- attr(Terms, "specials")$tt
if (missing(tt))
tt <- NULL
if (length(timetrans)) {
timetrans <- untangle.specials(Terms, "tt")
ntrans <- length(timetrans$terms)
if (is.null(tt)) {
tt <- function(x, time, riskset, weights) {
obrien <- function(x) {
r <- rank(x)
(r - 0.5)/(0.5 + length(r) - r)
}
unlist(tapply(x, riskset, obrien))
}
}
if (is.function(tt))
tt <- list(tt)
if (is.list(tt)) {
if (any(!sapply(tt, is.function)))
stop("The tt argument must contain function or list of functions")
if (length(tt) != ntrans) {
if (length(tt) == 1) {
temp <- vector("list", ntrans)
for (i in 1:ntrans) temp[[i]] <- tt[[1]]
tt <- temp
}
else stop("Wrong length for tt argument")
}
}
else stop("The tt argument must contain a function or list of functions")
if (ncol(Y) == 2) {
if (length(strats) == 0) {
sorted <- order(-Y[, 1], Y[, 2])
newstrat <- rep.int(0L, nrow(Y))
newstrat[1] <- 1L
}
else {
sorted <- order(strats, -Y[, 1], Y[, 2])
newstrat <- as.integer(c(1, 1 * (diff(strats[sorted]) !=
0)))
}
if (storage.mode(Y) != "double")
storage.mode(Y) <- "double"
counts <- .Call(Ccoxcount1, Y[sorted, ], as.integer(newstrat))
tindex <- sorted[counts$index]
}
else {
if (length(strats) == 0) {
sort.end <- order(-Y[, 2], Y[, 3])
sort.start <- order(-Y[, 1])
newstrat <- c(1L, rep(0, nrow(Y) - 1))
}
else {
sort.end <- order(strats, -Y[, 2], Y[, 3])
sort.start <- order(strats, -Y[, 1])
newstrat <- c(1L, as.integer(diff(strats[sort.end]) !=
0))
}
if (storage.mode(Y) != "double")
storage.mode(Y) <- "double"
counts <- .Call(Ccoxcount2, Y, as.integer(sort.start -
1L), as.integer(sort.end - 1L), as.integer(newstrat))
tindex <- counts$index
}
Y <- Surv(rep(counts$time, counts$nrisk), counts$status)
type <- "right"
mf <- mf[tindex, ]
strats <- rep(1:length(counts$nrisk), counts$nrisk)
weights <- model.weights(mf)
if (!is.null(weights) && any(!is.finite(weights)))
stop("weights must be finite")
tcall <- attr(Terms, "variables")[timetrans$terms + 2]
pvars <- attr(Terms, "predvars")
pmethod <- sub("makepredictcall.", "", as.vector(methods("makepredictcall")))
for (i in 1:ntrans) {
newtt <- (tt[[i]])(mf[[timetrans$var[i]]], Y[, 1],
strats, weights)
mf[[timetrans$var[i]]] <- newtt
nclass <- class(newtt)
if (any(nclass %in% pmethod)) {
dummy <- as.call(list(as.name(class(newtt)[1]),
tcall[[i]][[2]]))
ptemp <- makepredictcall(newtt, dummy)
pvars[[timetrans$terms[i] + 2]] <- ptemp
}
}
attr(Terms, "predvars") <- pvars
}
cluster <- attr(Terms, "specials")$cluster
if (length(cluster)) {
robust <- TRUE
tempc <- untangle.specials(Terms, "cluster", 1:10)
ord <- attr(Terms, "order")[tempc$terms]
if (any(ord > 1))
stop("Cluster can not be used in an interaction")
cluster <- strata(mf[, tempc$vars], shortlabel = TRUE)
dropterms <- tempc$terms
xlevels <- .getXlevels(Terms[-tempc$terms], mf)
}
else {
if (missing(robust))
robust <- FALSE
xlevels <- .getXlevels(Terms, mf)
dropterms <- NULL
}
contrast.arg <- NULL
attr(Terms, "intercept") <- 1
stemp <- untangle.specials(Terms, "strata", 1)
hasinteractions <- FALSE
if (length(stemp$vars) > 0) {
for (i in stemp$vars) {
if (any(attr(Terms, "order")[attr(Terms, "factors")[i,
] > 0] > 1))
hasinteractions <- TRUE
}
if (!hasinteractions)
dropterms <- c(dropterms, stemp$terms)
}
if (length(dropterms)) {
Terms2 <- Terms[-dropterms]
X <- model.matrix(Terms2, mf, constrasts = contrast.arg)
temp <- attr(X, "assign")
shift <- sort(dropterms)
temp <- temp + 1 * (shift[1] <= temp)
if (length(shift) == 2)
temp + 1 * (shift[2] <= temp)
attr(X, "assign") <- temp
}
else X <- model.matrix(Terms, mf, contrasts = contrast.arg)
if (!all(is.finite(X)))
stop("data contains an infinite predictor")
Xatt <- attributes(X)
if (hasinteractions)
adrop <- c(0, untangle.specials(Terms, "strata")$terms)
else adrop <- 0
xdrop <- Xatt$assign %in% adrop
X <- X[, !xdrop, drop = FALSE]
attr(X, "assign") <- Xatt$assign[!xdrop]
attr(X, "contrasts") <- Xatt$contrasts
offset <- model.offset(mf)
if (is.null(offset) | all(offset == 0))
offset <- rep(0, nrow(mf))
else if (any(!is.finite(offset) | !is.finite(exp(offset))))
stop("offsets must lead to a finite risk score")
weights <- model.weights(mf)
if (!is.null(weights) && any(!is.finite(weights)))
stop("weights must be finite")
assign <- attrassign(X, Terms)
contr.save <- attr(X, "contrasts")
if (missing(init))
init <- NULL
else {
if (length(init) != ncol(X))
stop("wrong length for init argument")
temp <- X %*% init - sum(colMeans(X) * init) + offset
if (any(exp(temp) > .Machine$double.xmax) || all(exp(temp) ==
0))
stop("initial values lead to overflow or underflow of the exp function")
}
if (sum(Y[, ncol(Y)]) == 0) {
ncoef <- ncol(X)
ctemp <- rep(NA, ncoef)
names(ctemp) <- colnames(X)
concordance = c(concordant = 0, discordant = 0, tied.x = 0,
tied.y = 0, tied.xy = 0, concordance = NA, std = NA,
timefix = FALSE)
rval <- list(coefficients = ctemp, var = matrix(0, ncoef,
ncoef), loglik = c(0, 0), score = 0, iter = 0, linear.predictors = offset,
residuals = rep(0, data.n), means = colMeans(X),
method = method, n = data.n, nevent = 0, terms = Terms,
assign = assign, concordance = concordance, y = Y,
call = Call)
class(rval) <- "coxph"
return(rval)
}
pterms <- sapply(mf, inherits, "coxph.penalty")
if (any(pterms)) {
pattr <- lapply(mf[pterms], attributes)
pname <- names(pterms)[pterms]
ord <- attr(Terms, "order")[match(pname, attr(Terms,
"term.labels"))]
if (any(ord > 1))
stop("Penalty terms cannot be in an interaction")
pcols <- assign[match(pname, names(assign))]
fit <- coxpenal.fit(X, Y, strats, offset, init = init,
control, weights = weights, method = method, row.names(mf),
pcols, pattr, assign)
}
else {
if (method == "breslow" || method == "efron") {
if (type == "right")
fitter <- get("coxph.fit")
else fitter <- get("agreg.fit")
}
else if (method == "exact") {
if (type == "right")
fitter <- get("coxexact.fit")
else fitter <- get("agexact.fit")
}
else stop(paste("Unknown method", method))
fit <- fitter(X, Y, strats, offset, init, control, weights = weights,
method = method, row.names(mf))
}
if (is.character(fit)) {
fit <- list(fail = fit)
class(fit) <- "coxph"
}
else {
if (!is.null(fit$coefficients) && any(is.na(fit$coefficients))) {
vars <- (1:length(fit$coefficients))[is.na(fit$coefficients)]
msg <- paste("X matrix deemed to be singular; variable",
paste(vars, collapse = " "))
if (!singular.ok)
stop(msg)
}
fit$n <- data.n
fit$nevent <- sum(Y[, ncol(Y)])
fit$terms <- Terms
fit$assign <- assign
class(fit) <- fit$class
fit$class <- NULL
if (robust && !is.null(fit$coefficients) && !all(is.na(fit$coefficients))) {
fit$naive.var <- fit$var
fit2 <- c(fit, list(x = X, y = Y, weights = weights))
if (length(strats))
fit2$strata <- strats
if (length(cluster)) {
temp <- residuals.coxph(fit2, type = "dfbeta",
collapse = cluster, weighted = TRUE)
if (is.null(init))
fit2$linear.predictors <- 0 * fit$linear.predictors
else fit2$linear.predictors <- c(X %*% init)
temp0 <- residuals.coxph(fit2, type = "score",
collapse = cluster, weighted = TRUE)
}
else {
temp <- residuals.coxph(fit2, type = "dfbeta",
weighted = TRUE)
fit2$linear.predictors <- 0 * fit$linear.predictors
temp0 <- residuals.coxph(fit2, type = "score",
weighted = TRUE)
}
fit$var <- t(temp) %*% temp
u <- apply(as.matrix(temp0), 2, sum)
fit$rscore <- coxph.wtest(t(temp0) %*% temp0, u,
control$toler.chol)$test
}
if (length(fit$coefficients) && is.null(fit$wald.test)) {
nabeta <- !is.na(fit$coefficients)
if (is.null(init))
temp <- fit$coefficients[nabeta]
else temp <- (fit$coefficients - init[1:length(fit$coefficients)])[nabeta]
fit$wald.test <- coxph.wtest(fit$var[nabeta, nabeta],
temp, control$toler.chol)$test
}
if (length(cluster))
temp <- concordancefit(Y, fit$linear.predictors,
strats, weights, cluster = cluster, reverse = TRUE,
timefix = FALSE)
else temp <- concordancefit(Y, fit$linear.predictors,
strats, weights, reverse = TRUE, timefix = FALSE)
if (is.matrix(temp$count))
fit$concordance <- c(colSums(temp$count), concordance = temp$concordance,
std = sqrt(temp$var))
else fit$concordance <- c(temp$count, concordance = temp$concordance,
std = sqrt(temp$var))
na.action <- attr(mf, "na.action")
if (length(na.action))
fit$na.action <- na.action
if (model) {
if (length(timetrans)) {
mf[[".surv."]] <- Y
mf[[".strata."]] <- strats
stop("'model=TRUE' not supported for models with tt terms")
}
fit$model <- mf
}
if (x) {
fit$x <- X
if (length(strats)) {
if (length(timetrans))
fit$strata <- strats
else fit$strata <- strata.keep
}
}
if (y)
fit$y <- Y
}
if (!is.null(weights) && any(weights != 1))
fit$weights <- weights
names(fit$means) <- names(fit$coefficients)
fit$formula <- formula(Terms)
if (length(xlevels) > 0)
fit$xlevels <- xlevels
fit$contrasts <- contr.save
if (any(offset != 0))
fit$offset <- offset
fit$call <- Call
fit
}
terms.inner<-function (x)
{
if (inherits(x, "formula")) {
if (length(x) == 3)
c(terms.inner(x[[2]]), terms.inner(x[[3]]))
else terms.inner(x[[2]])
}
else if (inherits(x, "call") && (x[[1]] != as.name("$") &&
x[[1]] != as.name("["))) {
if (x[[1]] == "+" || x[[1]] == "*" || x[[1]] == "-") {
if (length(x) == 3)
c(terms.inner(x[[2]]), terms.inner(x[[3]]))
else terms.inner(x[[2]])
}
else if (x[[1]] == as.name("Surv"))
unlist(lapply(x[-1], terms.inner))
else terms.inner(x[[2]])
}
else (deparse(x))
}
coxpenal.df<-function (hmat, hinv, fdiag, assign.list, ptype, nvar, pen1,
pen2, sparse)
{
if (ptype == 1 & nvar == 0) {
hdiag <- 1/fdiag
list(fvar2 = (hdiag - pen1) * fdiag^2, df = sum((hdiag -
pen1) * fdiag), fvar = fdiag, trH = sum(fdiag))
}
else if (ptype == 2) {
hmat.full <- t(hmat) %*% (ifelse(fdiag == 0, 0, 1/fdiag) *
hmat)
hinv.full <- hinv %*% (fdiag * t(hinv))
if (length(pen2) == length(hmat.full))
imat <- hmat.full - pen2
else imat <- hmat.full - diag(pen2)
var <- hinv.full %*% imat %*% hinv.full
if (length(assign.list) == 1)
list(var2 = var, df = sum(imat * hinv.full), trH = sum(diag(hinv.full)),
var = hinv.full)
else {
df <- trH <- NULL
d2 <- diag(hinv.full)
for (i in assign.list) {
temp <- coxph.wtest(hinv.full[i, i], var[i, i])$solve
if (is.matrix(temp))
df <- c(df, sum(diag(temp)))
else df <- c(df, sum(temp))
trH <- c(trH, sum(d2[i]))
}
list(var2 = var, df = df, trH = trH, var = hinv.full)
}
}
else {
nf <- length(fdiag) - nvar
nr1 <- 1:nf
nr2 <- (nf + 1):(nf + nvar)
d1 <- fdiag[nr1]
d2 <- fdiag[nr2]
temp <- t(hinv[nr1, ])
temp2 <- t(hinv[nr2, , drop = FALSE])
A.diag <- d1 + c(rep(1, nvar) %*% (temp^2 * d2))
B <- hinv[nr1, ] %*% (d2 * temp2)
C <- hinv[nr2, ] %*% (d2 * temp2)
var2 <- C - t(B) %*% (pen1 * B)
if (ptype == 3) {
hmat.22 <- t(hmat) %*% (ifelse(fdiag == 0, 0, 1/fdiag) *
hmat)
temp <- C - coxph.wtest(hmat.22, diag(nvar))$solve
if (nvar == 1) {
var2 <- var2 - C * pen2 * C
temp2 <- c(temp * pen2)
}
else if (length(pen2) == nvar) {
var2 <- var2 - C %*% (pen2 * C)
temp2 <- sum(diag(temp) * pen2)
}
else {
var2 <- var2 - C %*% matrix(pen2, nvar) %*% C
temp2 <- sum(diag(temp * pen2))
}
}
else temp2 <- 0
df <- trH <- NULL
cdiag <- diag(C)
for (i in 1:length(assign.list)) {
if (sparse == i) {
df <- c(df, nf - (sum(A.diag * pen1) + temp2))
trH <- c(trH, sum(A.diag))
}
else {
j <- assign.list[[i]]
temp <- coxph.wtest(C[j, j], var2[j, j])$solve
if (is.matrix(temp))
df <- c(df, sum(diag(temp)))
else df <- c(df, sum(temp))
trH <- c(trH, sum(cdiag[j]))
}
}
list(var = C, df = df, trH = trH, fvar = A.diag, var2 = var2)
}
}
coxpenal.fit<-function (x, y, strata, offset, init, control, weights, method,
rownames, pcols, pattr, assign)
{
eps <- control$eps
n <- nrow(y)
if (is.matrix(x))
nvar <- ncol(x)
else if (length(x) == 0)
stop("Must have an X variable")
else nvar <- 1
if (missing(offset) || is.null(offset))
offset <- rep(0, n)
if (missing(weights) || is.null(weights))
weights <- rep(1, n)
else {
if (any(weights <= 0))
stop("Invalid weights, must be >0")
}
if (ncol(y) == 3) {
if (length(strata) == 0) {
sorted <- cbind(order(-y[, 2], y[, 3]), order(-y[,
1])) - 1L
newstrat <- as.integer(n)
}
else {
sorted <- cbind(order(strata, -y[, 2], y[, 3]), order(strata,
-y[, 1])) - 1L
newstrat <- as.integer(cumsum(table(strata)))
}
status <- y[, 3]
andersen <- TRUE
}
else {
if (length(strata) == 0) {
sorted <- order(-y[, 1], y[, 2]) - 1L
newstrat <- as.integer(n)
}
else {
sorted <- order(strata, -y[, 1], y[, 2]) - 1L
newstrat <- as.integer(cumsum(table(strata)))
}
status <- y[, 2]
andersen <- FALSE
}
n.eff <- sum(y[, ncol(y)])
npenal <- length(pattr)
if (npenal == 0 || length(pcols) != npenal)
stop("Invalid pcols or pattr arg")
sparse <- sapply(pattr, function(x) !is.null(x$sparse) &&
x$sparse)
if (sum(sparse) > 1)
stop("Only one sparse penalty term allowed")
pterms <- rep(0, length(assign))
names(pterms) <- names(assign)
pindex <- rep(0, npenal)
for (i in 1:npenal) {
temp <- unlist(lapply(assign, function(x, y) (length(x) ==
length(y) && all(x == y)), pcols[[i]]))
if (sparse[i])
pterms[temp] <- 2
else pterms[temp] <- 1
pindex[i] <- (seq(along.with = temp))[temp]
}
if ((sum(pterms == 2) != sum(sparse)) || (sum(pterms > 0) !=
npenal))
stop("pcols and assign arguments disagree")
if (any(pindex != sort(pindex))) {
temp <- order(pindex)
pindex <- pindex[temp]
pcols <- pcols[temp]
pattr <- pattr[temp]
}
ptype <- any(sparse) + 2 * (any(!sparse))
f.expr1 <- function(coef) NULL
f.expr2 <- function(coef) NULL
if (any(sparse)) {
sparse.attr <- (pattr[sparse])[[1]]
fcol <- unlist(pcols[sparse])
if (length(fcol) > 1)
stop("Sparse term must be single column")
xx <- x[, -fcol, drop = FALSE]
for (i in 1:length(assign)) {
j <- assign[[i]]
if (j[1] > fcol)
assign[[i]] <- j - 1
}
for (i in 1:npenal) {
j <- pcols[[i]]
if (j[1] > fcol)
pcols[[i]] <- j - 1
}
frailx <- x[, fcol]
frailx <- match(frailx, sort(unique(frailx)))
nfrail <- max(frailx)
nvar <- nvar - 1
pfun1 <- sparse.attr$pfun
f.expr1 <- function(coef) {
coxlist1$coef <- coef
if (is.null(extra1))
temp <- pfun1(coef, theta1, n.eff)
else temp <- pfun1(coef, theta1, n.eff, extra1)
if (!is.null(temp$recenter))
coxlist1$coef <- coxlist1$coef - as.double(temp$recenter)
if (!temp$flag) {
coxlist1$first <- -as.double(temp$first)
coxlist1$second <- as.double(temp$second)
}
coxlist1$penalty <- -as.double(temp$penalty)
coxlist1$flag <- as.logical(temp$flag)
if (any(sapply(coxlist1, length) != c(rep(nfrail,
3), 1, 1)))
stop("Incorrect length in coxlist1")
coxlist1
}
if (!is.null(getOption("survdebug")))
debug(f.expr1)
coxlist1 <- list(coef = double(nfrail), first = double(nfrail),
second = double(nfrail), penalty = 0, flag = FALSE)
}
else {
xx <- x
frailx <- 0
nfrail <- 0
}
if (sum(!sparse) > 0) {
full.imat <- !all(unlist(lapply(pattr, function(x) x$diag)))
ipenal <- (1:length(pattr))[!sparse]
f.expr2 <- function(coef) {
coxlist2$coef <- coef
pentot <- 0
for (i in ipenal) {
pen.col <- pcols[[i]]
coef <- coxlist2$coef[pen.col]
if (is.null(extralist[[i]]))
temp <- ((pattr[[i]])$pfun)(coef, thetalist[[i]],
n.eff)
else temp <- ((pattr[[i]])$pfun)(coef, thetalist[[i]],
n.eff, extralist[[i]])
if (!is.null(temp$recenter))
coxlist2$coef[pen.col] <- coxlist2$coef[pen.col] -
temp$recenter
if (temp$flag)
coxlist2$flag[pen.col] <- TRUE
else {
coxlist2$flag[pen.col] <- FALSE
coxlist2$first[pen.col] <- -temp$first
if (full.imat) {
tmat <- matrix(coxlist2$second, nvar, nvar)
tmat[pen.col, pen.col] <- temp$second
coxlist2$second <- c(tmat)
}
else coxlist2$second[pen.col] <- temp$second
}
pentot <- pentot - temp$penalty
}
coxlist2$penalty <- as.double(pentot)
if (any(sapply(coxlist2, length) != length2))
stop("Length error in coxlist2")
coxlist2
}
if (!is.null(getOption("survdebug")))
debug(f.expr2)
if (full.imat) {
coxlist2 <- list(coef = double(nvar), first = double(nvar),
second = double(nvar * nvar), penalty = 0, flag = rep(FALSE,
nvar))
length2 <- c(nvar, nvar, nvar * nvar, 1, nvar)
}
else {
coxlist2 <- list(coef = double(nvar), first = double(nvar),
second = double(nvar), penalty = 0, flag = rep(FALSE,
nvar))
length2 <- c(nvar, nvar, nvar, 1, nvar)
}
}
else full.imat <- FALSE
if (nfrail > 0)
finit <- rep(0, nfrail)
else finit <- 0
if (!missing(init) && !is.null(init)) {
if (length(init) != nvar) {
if (length(init) == (nvar + nfrail)) {
finit <- init[-(1:nvar)]
init <- init[1:nvar]
}
else stop("Wrong length for inital values")
}
}
else init <- double(nvar)
cfun <- lapply(pattr, function(x) x$cfun)
parmlist <- lapply(pattr, function(x, eps) c(x$cparm, eps2 = eps),
sqrt(eps))
extralist <- lapply(pattr, function(x) x$pparm)
iterlist <- vector("list", length(cfun))
thetalist <- vector("list", length(cfun))
printfun <- lapply(pattr, function(x) x$printfun)
for (i in 1:length(cfun)) {
temp <- (cfun[[i]])(parmlist[[i]], iter = 0)
if (sparse[i]) {
theta1 <- temp$theta
extra1 <- extralist[[i]]
}
thetalist[[i]] <- temp$theta
iterlist[[i]] <- temp
}
temp1 <- c("x", "coef", "plik", "loglik", "status", "neff",
"df", "trH")
temp2 <- c("frailx", "coxfit$fcoef", "loglik1", "coxfit$loglik",
"status", "n.eff")
temp3 <- c("xx[,pen.col]", "coxfit$coef[pen.col]", "loglik1",
"coxfit$loglik", "status", "n.eff")
calls <- vector("expression", length(cfun))
cargs <- lapply(pattr, function(x) x$cargs)
for (i in 1:length(cfun)) {
tempchar <- paste("(cfun[[", i, "]])(parmlist[[", i,
"]], iter,", "iterlist[[", i, "]]")
temp2b <- c(temp2, paste("pdf[", i, "]"), paste("trH[",
i, "]"))
temp3b <- c(temp3, paste("pdf[", i, "]"), paste("trH[",
i, "]"))
if (length(cargs[[i]]) == 0)
calls[i] <- parse(text = paste(tempchar, ")"))
else {
temp <- match(cargs[[i]], temp1)
if (any(is.na(temp)))
stop(paste((cargs[[i]])[is.na(temp)], "not matched"))
if (sparse[i])
temp4 <- paste(temp2b[temp], collapse = ",")
else temp4 <- paste(temp3b[temp], collapse = ",")
calls[i] <- parse(text = paste(paste(tempchar, temp4,
sep = ","), ")"))
}
}
need.df <- any(!is.na(match(c("df", "trH"), unlist(cargs))))
varnames <- dimnames(xx)[[2]]
for (i in 1:length(cfun)) {
if (!is.null(pattr[[i]]$varname))
varnames[pcols[[i]]] <- pattr[[i]]$varname
}
rho <- environment()
storage.mode(y) <- storage.mode(weights) <- "double"
storage.mode(xx) <- storage.mode(offset) <- "double"
if (andersen)
coxfit <- .C(Cagfit5a, as.integer(n), as.integer(nvar),
y, xx, offset, weights, newstrat, sorted, means = double(nvar),
coef = as.double(init), u = double(nvar), loglik = double(1),
as.integer(method == "efron"), as.integer(ptype),
as.integer(full.imat), as.integer(nfrail), as.integer(frailx),
f.expr1, f.expr2, rho)
else coxfit <- .C(Ccoxfit5a, as.integer(n), as.integer(nvar),
y, xx, offset, weights, newstrat, sorted, means = double(nvar),
coef = as.double(init), u = double(nvar), loglik = double(1),
as.integer(method == "efron"), as.integer(ptype), as.integer(full.imat),
as.integer(nfrail), as.integer(frailx), f.expr1, f.expr2,
rho)
loglik0 <- coxfit$loglik
means <- coxfit$means
iter2 <- 0
iterfail <- NULL
thetasave <- unlist(thetalist)
for (outer in 1:control$outer.max) {
if (andersen)
coxfit <- .C(Cagfit5b, iter = as.integer(control$iter.max),
as.integer(n), as.integer(nvar), as.integer(newstrat),
coef = as.double(init), u = double(nvar + nfrail),
hmat = double(nvar * (nvar + nfrail)), hinv = double(nvar *
(nvar + nfrail)), loglik = double(1), flag = integer(1),
as.double(control$eps), as.double(control$toler.chol),
as.integer(method == "efron"), as.integer(nfrail),
fcoef = as.double(finit), fdiag = double(nfrail +
nvar), f.expr1, f.expr2, rho)
else coxfit <- .C(Ccoxfit5b, iter = as.integer(control$iter.max),
as.integer(n), as.integer(nvar), as.integer(newstrat),
coef = as.double(init), u = double(nvar + nfrail),
hmat = double(nvar * (nvar + nfrail)), hinv = double(nvar *
(nvar + nfrail)), loglik = double(1), flag = integer(1),
as.double(control$eps), as.double(control$toler.chol),
as.integer(method == "efron"), as.integer(nfrail),
fcoef = as.double(finit), fdiag = double(nfrail +
nvar), f.expr1, f.expr2, rho)
iter <- outer
iter2 <- iter2 + coxfit$iter
if (coxfit$iter >= control$iter.max)
iterfail <- c(iterfail, iter)
temp <- rep(FALSE, nvar + nfrail)
if (nfrail > 0)
temp[1:nfrail] <- coxlist1$flag
if (ptype > 1)
temp[nfrail + 1:nvar] <- coxlist2$flag
fdiag <- ifelse(temp, 0, coxfit$fdiag)
if (need.df) {
if (nfrail > 0)
temp1 <- coxlist1$second
else temp1 <- 0
if (ptype > 1)
temp2 <- coxlist2$second
else temp2 <- 0
dftemp <- coxpenal.df(matrix(coxfit$hmat, ncol = nvar),
matrix(coxfit$hinv, ncol = nvar), fdiag, assign,
ptype, nvar, temp1, temp2, pindex[sparse])
df <- dftemp$df
var <- dftemp$var
var2 <- dftemp$var2
pdf <- df[pterms > 0]
trH <- dftemp$trH[pterms > 0]
}
if (nfrail > 0)
penalty <- -coxlist1$penalty
else penalty <- 0
if (ptype > 1)
penalty <- penalty - coxlist2$penalty
loglik1 <- coxfit$loglik + penalty
if (iter == 1)
penalty0 <- penalty
done <- TRUE
for (i in 1:length(cfun)) {
pen.col <- pcols[[i]]
temp <- eval(calls[i])
if (sparse[i])
theta1 <- temp$theta
thetalist[[i]] <- temp$theta
iterlist[[i]] <- temp
done <- done & temp$done
}
if (done)
break
if (iter == 1) {
init <- coefsave <- coxfit$coef
finit <- fsave <- coxfit$fcoef
thetasave <- cbind(thetasave, unlist(thetalist))
}
else {
temp <- as.vector(unlist(thetalist))
coefsave <- cbind(coefsave, coxfit$coef)
fsave <- cbind(fsave, coxfit$fcoef)
howclose <- apply((thetasave - temp)^2, 2, sum)
which <- min((1:iter)[howclose == min(howclose)])
if (nvar > 0)
init <- coefsave[, which]
if (nfrail > 0)
finit <- fsave[, which]
thetasave <- cbind(thetasave, temp)
}
}
if (nfrail > 0) {
lp <- offset + coxfit$fcoef[frailx]
if (nvar > 0)
lp <- lp + x[, -fcol, drop = FALSE] %*% coxfit$coef -
sum(means * coxfit$coef)
}
else lp <- offset + as.vector(x %*% coxfit$coef) - sum(means *
coxfit$coef)
if (andersen) {
.C(Cagfit5c, as.integer(nvar))
resid <- .Call(Cagmart3, y, exp(lp), weights, newstrat,
sorted, as.integer(method == "efron"))
}
else {
expect <- .C(Ccoxfit5c, as.integer(n), as.integer(nvar),
as.integer(newstrat), as.integer(method == "efron"),
expect = double(n))$expect
resid <- status - expect
}
names(resid) <- rownames
if (!need.df) {
if (nfrail > 0)
temp1 <- coxlist1$second
else temp1 <- 0
if (ptype > 1)
temp2 <- coxlist2$second
else temp2 <- 0
dftemp <- coxpenal.df(matrix(coxfit$hmat, ncol = nvar),
matrix(coxfit$hinv, ncol = nvar), fdiag, assign,
ptype, nvar, temp1, temp2, pindex[sparse])
df <- dftemp$df
trH <- dftemp$trH
var <- dftemp$var
var2 <- dftemp$var2
}
if (control$iter.max > 1 && length(iterfail) > 0)
warning(paste("Inner loop failed to coverge for iterations",
paste(iterfail, collapse = " ")))
which.sing <- (fdiag[nfrail + 1:nvar] == 0)
coef <- coxfit$coef
names(coef) <- varnames
coef[which.sing] <- NA
names(iterlist) <- names(pterms[pterms > 0])
if (nfrail > 0) {
if (nvar > 0) {
list(coefficients = coef, var = var, var2 = var2,
loglik = c(loglik0, loglik1), iter = c(iter,
iter2), linear.predictors = as.vector(lp),
residuals = resid, means = means, method = method,
class = c("coxph.penal", "coxph"), frail = coxfit$fcoef,
fvar = dftemp$fvar, df = df, df2 = dftemp$df2,
penalty = c(penalty0, penalty), pterms = pterms,
assign2 = assign, history = iterlist, coxlist1 = coxlist1,
printfun = printfun)
}
else {
list(loglik = c(loglik0, loglik1), iter = c(iter,
iter2), linear.predictors = as.vector(lp), residuals = resid,
means = means, method = method, class = c("coxph.penal",
"coxph"), frail = coxfit$fcoef, fvar = dftemp$fvar,
df = df, df2 = dftemp$df2, penalty = c(penalty0,
penalty), pterms = pterms, assign2 = assign,
history = iterlist, printfun = printfun)
}
}
else {
list(coefficients = coef, var = var, var2 = var2, loglik = c(loglik0,
loglik1), iter = c(iter, iter2), linear.predictors = lp,
residuals = resid, means = means, method = method,
class = c("coxph.penal", "coxph"), df = df, df2 = dftemp$df2,
penalty = c(penalty0, penalty), pterms = pterms,
assign2 = assign, history = iterlist, coxlist2 = coxlist2,
printfun = printfun)
}
}
parsecovar1<-function (flist, statedatanames)
{
if (any(sapply(flist, function(x) !inherits(x, "formula"))))
stop("flist must be a list of formulas")
if (any(sapply(flist, length) != 3))
stop("all formulas must have a left and right side")
lhs <- lapply(flist, function(x) x[-3])
rhs <- lapply(flist, function(x) x[-2])
rh2 <- lapply(rhs, function(form) {
parts <- strsplit(deparse(form, width.cutoff = 300, control = NULL),
"/", fixed = TRUE)[[1]]
if (length(parts) == 1) {
ival <- NULL
common <- FALSE
fixed <- FALSE
clear <- FALSE
}
else {
optterms <- terms(formula(paste("~", parts[2])))
ff <- rownames(attr(optterms, "factors"))
index <- match(ff, c("common", "fixed", "init", "clear"))
if (any(is.na(index)))
stop("option not recognized in a covariates formula: ",
paste(ff[is.na(index)], collapse = ", "))
common <- any(index == 1)
fixed <- any(index == 2)
clear <- any(index == 3)
if (any(index == 3)) {
optatt <- attributes(optterms)
j <- optatt$variables[1 + which(index == 3)]
j[[1]] <- as.name("list")
ival <- unlist(eval(j, parent.frame()))
}
else ival <- NULL
}
form <- formula(paste("~ -1 +", parts[1]))
list(common = common, fixed = fixed, clear = clear, ival = ival,
formula = form)
})
pcut <- function(form) {
if (length(form) == 3) {
if (form[[1]] == "+")
c(pcut(form[[2]]), pcut(form[[3]]))
else if (form[[1]] == "~")
pcut(form[[2]])
else list(form)
}
else list(form)
}
lcut <- lapply(lhs, function(x) pcut(x[[2]]))
env1 <- new.env(parent = parent.frame(2))
env2 <- new.env(parent = env1)
if (missing(statedatanames)) {
assign("state", function(...) list(stateid = "state",
values = c(...)), env1)
assign("state", list(stateid = "state"))
}
else {
for (i in statedatanames) {
assign(i, eval(list(stateid = i)), env2)
tfun <- eval(parse(text = paste0("function(...) list(stateid='",
i, "', values=c(...)")))
assign(i, tfun, env1)
}
}
lterm <- lapply(lcut, function(x) {
lapply(x, function(z) {
if (length(z) == 1) {
temp <- eval(z, envir = env2)
if (is.list(temp) && names(temp)[[1]] == "stateid")
temp
else temp
}
else if (length(z) == 3 && z[[1]] == ":")
list(left = eval(z[[2]], envir = env2), right = eval(z[[3]],
envir = env2))
else stop("invalid term: ", deparse(z))
})
})
list(rhs = rh2, lhs = lterm)
}
parsecovar2<-function (covar1, statedata, dformula, Terms, transitions, states)
{
if (is.null(statedata))
statedata <- data.frame(state = states)
else {
if (is.null(statedata$state))
stop("the statedata data set must contain a variable 'state'")
indx1 <- match(states, statedata$state, nomatch = 0)
if (any(indx1 == 0))
stop("statedata does not contain all the possible states: ",
states[indx1 == 0])
statedata <- statedata[indx1, ]
}
allterm <- attr(Terms, "term.labels")
nterm <- length(allterm)
nstate <- length(states)
tmap <- array(0L, dim = c(nterm + 1, nstate, nstate))
dterms <- match(attr(terms.formula(dformula), "term.labels"),
allterm)
dterms <- c(1L, 1L + dterms)
k <- seq(along = dterms)
for (i in 1:nstate) {
for (j in 1:nstate) {
tmap[dterms, i, j] <- k
k <- k + length(k)
}
}
ncoef <- max(tmap)
inits <- NULL
if (!is.null(covar1)) {
for (i in 1:length(covar1$rhs)) {
rhs <- covar1$rhs[[i]]
lhs <- covar1$lhs[[i]]
rterm <- terms.formula(rhs$formula)
rindex <- 1L + match(attr(rterm, "term.labels"),
allterm, nomatch = 0)
if (any(rindex == 1L))
stop("dterm mismatch bug 2")
if (attr(rterm, "intercept") == 1)
rindex <- c(1L, rindex)
state1 <- state2 <- NULL
for (x in lhs) {
if (is.null(x$left))
stop("term found without a :", x)
if (!is.list(x$left) && length(x$left) == 1 &
x$left == 1)
temp1 <- 1:nrow(statedata)
else if (is.list(x$left) && names(x$left)[1] ==
"stateid") {
if (is.null(x$left$value))
stop("state variable with no list of values: ",
x$left$stateid)
else temp1 <- which(statedata[[x$left$stateid]] %in%
x$left$value)
}
else temp1 <- which(statedata$state %in% x$left)
if (!is.list(x$right) && length(x$right) == 1 &&
x$right == 1)
temp2 <- 1:nrow(statedata)
else if (is.list(x$right) && names(x$right)[1] ==
"stateid") {
if (is.null(x$right$value))
stop("state variable with no list of values: ",
x$right$stateid)
else temp2 <- which(statedata[[x$right$stateid]] %in%
x$right$value)
}
else temp2 <- which(statedata$state %in% x$right)
state1 <- c(state1, rep(temp1, length(temp2)))
state2 <- c(state2, rep(temp2, each = length(temp1)))
}
if (rhs$clear)
tmap[-1, state1, state2] <- 0
if (length(rhs$ival))
inits <- c(inits, list(term = rindex, state1 = state1,
state2 = state2, init = rhs$ival))
if (rhs$common)
j <- ncoef + seq_len(length(rindex))
else j <- ncoef + seq_len(length(rindex) * length(state1))
tmap[rindex, state1, state2] <- j
ncoef <- max(j)
}
}
t2 <- transitions[, -1, drop = FALSE]
indx1 <- match(rownames(t2), states)
indx2 <- match(colnames(t2), states)
tmap2 <- matrix(0L, nrow = 1 + nterm, ncol = sum(t2 > 0))
trow <- row(t2)[t2 > 0]
tcol <- col(t2)[t2 > 0]
for (i in 1:length(trow)) tmap2[, i] <- tmap[, indx1[trow[i]],
indx2[tcol[i]]]
dimnames(tmap2) <- list(c("Intercept", allterm), paste(indx1[trow],
indx2[tcol], sep = ":"))
list(tmap = tmap2, inits = inits, mapid = cbind(indx1[trow],
indx2[tcol]))
}
parsecovar3<-function (tmap, Xcol, Xassign)
{
hasintercept <- (Xassign[1] == 0)
cmap <- matrix(0L, length(Xcol) + !hasintercept, ncol(tmap))
cmap[1, ] <- match(tmap[1, ], sort(c(0, unique(tmap[1, ])))) -
1L
xcount <- table(factor(Xassign, levels = 1:max(Xassign)))
mult <- 1 + max(xcount)
j <- 1
for (i in 2:nrow(tmap)) {
k <- seq_len(xcount[i - 1])
cmap[j + k, ] <- ifelse(tmap[i, ] == 0, 0, tmap[i, ] *
mult + rep(k, ncol(tmap)))
j <- j + max(k)
}
cmap[-1, ] <- match(cmap[-1, ], sort(unique(c(0L, cmap[-1,
])))) - 1L
colnames(cmap) <- colnames(tmap)
if (hasintercept)
rownames(cmap) <- Xcol
else rownames(cmap) <- c("(Intercept)", Xcol)
cmap
}
residuals.coxph<-function (object, type = c("martingale", "deviance", "score",
"schoenfeld", "dfbeta", "dfbetas", "scaledsch", "partial"),
collapse = FALSE, weighted = FALSE, ...)
{
type <- match.arg(type)
otype <- type
if (type == "dfbeta" || type == "dfbetas") {
otype <- type
type <- "score"
if (missing(weighted))
weighted <- TRUE
}
if (type == "scaledsch")
type <- "schoenfeld"
n <- length(object$residuals)
rr <- object$residuals
y <- object$y
x <- object[["x"]]
vv <- drop(object$naive.var)
if (is.null(vv))
vv <- drop(object$var)
weights <- object$weights
if (is.null(weights))
weights <- rep(1, n)
strat <- object$strata
method <- object$method
if (method == "exact" && (type == "score" || type == "schoenfeld"))
stop(paste(otype, "residuals are not available for the exact method"))
if (type == "martingale" || type == "partial")
rr <- object$residuals
else {
Terms <- object$terms
if (!inherits(Terms, "terms"))
stop("invalid terms component of object")
strats <- attr(Terms, "specials")$strata
if (is.null(y) || (is.null(x) && type != "deviance")) {
temp <- coxph.getdata(object, y = TRUE, x = TRUE,
stratax = TRUE)
y <- temp$y
x <- temp$x
if (length(strats))
strat <- temp$strata
}
ny <- ncol(y)
status <- y[, ny, drop = TRUE]
if (type != "deviance") {
nstrat <- as.numeric(strat)
nvar <- ncol(x)
if (is.null(strat)) {
ord <- order(y[, ny - 1], -status)
newstrat <- rep(0, n)
}
else {
ord <- order(nstrat, y[, ny - 1], -status)
newstrat <- c(diff(as.numeric(nstrat[ord])) !=
0, 1)
}
newstrat[n] <- 1
x <- x[ord, ]
y <- y[ord, ]
score <- exp(object$linear.predictors)[ord]
}
}
if (type == "schoenfeld") {
if (ny == 2) {
mintime <- min(y[, 1])
if (mintime < 0)
y <- cbind(2 * mintime - 1, y)
else y <- cbind(-1, y)
}
temp <- .C(Ccoxscho, n = as.integer(n), as.integer(nvar),
as.double(y), resid = as.double(x), as.double(score *
weights[ord]), as.integer(newstrat), as.integer(method ==
"efron"), double(3 * nvar))
deaths <- y[, 3] == 1
if (nvar == 1)
rr <- temp$resid[deaths]
else rr <- matrix(temp$resid[deaths], ncol = nvar)
if (weighted)
rr <- rr * weights[deaths]
if (length(strats))
attr(rr, "strata") <- table((strat[ord])[deaths])
time <- c(y[deaths, 2])
if (is.matrix(rr))
dimnames(rr) <- list(time, names(object$coefficients))
else names(rr) <- time
if (otype == "scaledsch") {
ndead <- sum(deaths)
coef <- ifelse(is.na(object$coefficients), 0, object$coefficients)
rr <- drop(rr %*% vv * ndead + rep(coef, each = nrow(rr)))
}
return(rr)
}
if (type == "score") {
if (ny == 2) {
resid <- .C(Ccoxscore, as.integer(n), as.integer(nvar),
as.double(y), x = as.double(x), as.integer(newstrat),
as.double(score), as.double(weights[ord]), as.integer(method ==
"efron"), resid = double(n * nvar), double(2 *
nvar))$resid
}
else {
resid <- .C(Cagscore, as.integer(n), as.integer(nvar),
as.double(y), as.double(x), as.integer(newstrat),
as.double(score), as.double(weights[ord]), as.integer(method ==
"efron"), resid = double(n * nvar), double(nvar *
6))$resid
}
if (nvar > 1) {
rr <- matrix(0, n, nvar)
rr[ord, ] <- matrix(resid, ncol = nvar)
dimnames(rr) <- list(names(object$residuals), names(object$coefficients))
}
else rr[ord] <- resid
if (otype == "dfbeta") {
if (is.matrix(rr))
rr <- rr %*% vv
else rr <- rr * vv
}
else if (otype == "dfbetas") {
if (is.matrix(rr))
rr <- (rr %*% vv) %*% diag(sqrt(1/diag(vv)))
else rr <- rr * sqrt(vv)
}
}
if (weighted)
rr <- rr * weights
if (!is.null(object$na.action)) {
rr <- naresid(object$na.action, rr)
if (is.matrix(rr))
n <- nrow(rr)
else n <- length(rr)
if (type == "deviance")
status <- naresid(object$na.action, status)
}
if (type == "partial") {
rr <- rr + predict(object, type = "terms")
}
if (!missing(collapse)) {
if (length(collapse) != n)
stop("Wrong length for 'collapse'")
rr <- drop(rowsum(rr, collapse))
if (type == "deviance")
status <- drop(rowsum(status, collapse))
}
if (type == "deviance")
sign(rr) * sqrt(-2 * (rr + ifelse(status == 0, 0, status *
log(status - rr))))
else rr
}
coxph.getdata<-function (fit, y = TRUE, x = TRUE, stratax = TRUE, offset = FALSE)
{
ty <- fit[["y"]]
tx <- fit[["x"]]
strat <- fit$strata
Terms <- fit$terms
if (is.null(attr(Terms, "offset")))
offset <- FALSE
if (offset)
x <- TRUE
if (!inherits(Terms, "terms"))
stop("invalid terms component of fit")
strats <- attr(Terms, "specials")$strata
if (length(strats) == 0)
stratax <- FALSE
if ((y && is.null(ty)) || (x && is.null(tx)) || (stratax &&
is.null(strat)) || offset) {
m <- stats::model.frame(fit)
if (y && is.null(ty))
ty <- model.extract(m, "response")
if (offset)
toff <- model.extract(m, "offset")
if ((x || stratax) && is.null(tx)) {
if (stratax) {
temp <- untangle.specials(Terms, "strata", 1)
strat <- strata(m[temp$vars], shortlabel = TRUE)
}
if (x)
tx <- model.matrix(fit, data = m)
}
}
else if (offset)
toff <- fit$linear.predictors - (c(tx %*% fit$coef) -
sum(fit$means * fit$coef))
temp <- list()
if (y)
temp$y <- ty
if (x)
temp$x <- tx
if (stratax)
temp$strata <- strat
if (offset)
temp$offset <- toff
temp
}
# Create a strata variable, possibly from many objects
#
strata <- function(..., na.group=FALSE, shortlabel, sep=', ') {
# First, grab a copy of the call, which will be used to manufacture
# labels for unlabeled arguments
# Then get the arguments as a list
words <- as.character((match.call())[-1])
allf <- list(...)
# If there is only one argument, and it itself is a list, use
# it instead
if(length(allf) == 1 && is.list(ttt <- unclass(allf[[1]]))) allf <- ttt
nterms <- length(allf)
# Keep the names of named args as their label, what was typed otherwise
if (is.null(names(allf))) {
argname <- words[1:nterms]
if (missing(shortlabel))
shortlabel <- all(sapply(allf,
function(x) is.character(x) | inherits(x, 'factor')))
}
else {
argname <- ifelse(names(allf) == '', words[1:nterms], names(allf))
if (missing(shortlabel)) shortlabel <- FALSE
}
# If the arguments are not all the same length, stop now
# Mostly this is to stop calls with an improper object
arglength <- sapply(allf, length)
if (any(arglength != arglength[1]))
stop("all arguments must be the same length")
if (!all(sapply(allf, is.atomic))) stop("all arguments must be vectors")
# Process the first argument
what <- allf[[1]]
if(is.null(levels(what)))
what <- factor(what)
levs <- unclass(what) - 1
wlab <- levels(what)
if (na.group && any(is.na(what))){
# add "NA" as a level
levs[is.na(levs)] <- length(wlab)
wlab <- c(wlab, "NA")
}
if (shortlabel) labs <- wlab
else labs <- paste(argname[1], wlab, sep='=')
# Now march through the other variables, if any
for (i in (1:nterms)[-1]) {
what <- allf[[i]]
if(is.null(levels(what)))
what <- factor(what)
wlab <- levels(what)
wlev <- unclass(what) - 1
if (na.group && any(is.na(wlev))){
wlev[is.na(wlev)] <- length(wlab)
wlab <- c(wlab, "NA")
}
if (!shortlabel) wlab <- format(paste(argname[i], wlab, sep='='))
levs <- wlev + levs*(length(wlab))
labs <- paste(rep(labs, rep(length(wlab), length(labs))),
rep(wlab, length(labs)), sep=sep)
}
levs <- levs + 1
ulevs <- sort(unique(levs[!is.na(levs)]))
levs <- match(levs, ulevs)
labs <- labs[ulevs]
factor(levs, labels=labs)
}
gerstung-lab/ebmstate documentation built on Feb. 11, 2022, 4:01 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions strata MakeInteger print.msfit print.coxrfx summary.coxrfx CoxRFX
Documented in CoxRFX MakeInteger print.coxrfx print.msfit summary.coxrfx
#' Empirical Bayes, multi-state Cox model
#'
#' This function estimates a multi-state Cox model with one or more Gaussian priors
#' imposed on the regression coefficients (see Therneau et al., 2003).
#' Multiple groups of coefficients can be defined: coefficients within a group share
#' the same (possibly unknown) mean and variance. The parameters and hyperparameters are
#' efficiently estimated by an EM-type algorithm built around the function
#' \code{survival::coxph}.
#' @param Z A data frame consisting of the covariate columns of a data set in 'long format',
#' and two extra columns: one named 'trans', with
#' the transition that each row refers to, and
#' another named 'strata', with the stratum
#' of each transition (transitions belonging to the same
#' stratum are assumed to have the same baseline hazard
#' function).
#' @param surv A `survival' object created with \code{survival::Surv}.
#' @param groups A character or numeric vector whose \eqn{i}th element gives the group of the regression
#' coefficient associated with the \eqn{i}th covariate column of Z (coefficients belonging to the same group
#' share the same Gaussian prior).
#' @param tmat Transition matrix describing the states and transitions in
#' the multi-state model. See \code{trans} in \code{\link[mstate:msprep]{mstate::msprep}}
#' for more detailed information.
#' @param which.mu A vector with names or numbers of coefficient groups (see
#' argument \code{groups}). If the name or number of a group of coefficients is
#' given in this argument, \code{CoxRFX} will estimate the mean of its Gaussian distribution;
#' otherwise the mean will be fixed at zero.
#' @param tol Convergence criterium of the EM algorithm. The algorithm stops unless
#' there is at least one parameter (or hyperparameter) for which it holds that the
#' current estimate differs in absolute terms by more than \code{tol} from the
#' previous estimate.
#' @param max.iter The maximum number of iterations in the EM algorithm.
#' @param sigma0 A vector with the initial value of the variance hyperparameter for each group of coefficients.
#' Or a single value, in case the initial value of the variance hyperparameter is meant to be the same for all groups.
#' @param sigma.hat Which estimator to use for the variance hyperparameters (see details).
#' @param verbose Gives more output.
#' @param ... Further arguments passed to the function \code{survival::coxph}.
#' @details
#'
#' Different estimators exist for the variance hyperparameters: the default is "df", as used by Perperoglou (2014) and introduced by Schall (1991).
#' Alternatives are MLE, REML, and BLUP, as defined by Therneau et al. (2003).
#' Simulations suggest that the 'df' method is the most accurate.
#'
#' The model can also be fitted using package \code{coxme}; the \code{coxme}
#' routine numerically optimises the integrated partial likelihood, which may
#' be more accurate, but is computationally expensive.
#'
#' @references Terry M Therneau, Patricia M Grambsch & V. Shane Pankratz (2003) Penalized Survival Models and Frailty, Journal of Computational and Graphical Statistics, 12:1, 156-175, http://dx.doi.org/10.1198/1061860031365
#'
#' A. Perperoglou (2014). Cox models with dynamic ridge penalties on time-varying effects of the covariates. Stat Med, 33:170-80. http://dx.doi.org/10.1002/sim.5921
#'
#' R. Schall (1991). Estimation in generalized linear models with random effects. Biometrika, 78:719-727. http://dx.doi.org/10.1093/biomet/78.4.719
#' @return An object of class \code{c(coxrfx,coxph.penal,coxph)},
#' which is essentially a \code{coxph} object with a few extra
#' fields [the inputs $groups, $Z, $surv, and $tmat,
#' and the hyperparameters $sigma2 (variances) and $mu (means)].
#' See \code{\link[survival:coxph.object]{survival::coxph.object}}.
#'
#' @author Moritz Gerstung & Rui Costa, extending the work of
#' Terry Therneau et al. in the package \code{survival}.
#' @seealso Package \code{survival}
#' \code{\link[survival:coxph.object]{survival::coxph.object}};
#' \code{\link[survival:Surv]{survival::Surv}}; package \code{coxme}.
#' @export
#' @example inst/examples/CoxRFX-example.R
CoxRFX <- function(Z, surv, groups = rep(1, ncol(Z)),tmat=NULL, which.mu = unique(groups), tol=1e-3, max.iter=50, sigma0 = 0.1, sigma.hat=c("df","MLE","REML","BLUP"), verbose=FALSE, ...){
##
trans<-Z$trans
strata<-Z$strata
Z$strata<-Z$trans<-NULL
namesZ<-names(Z)
##
Z = as.matrix(Z)
Z.df <- TRUE
if(is.null(colnames(Z)))
colnames(Z) <- make.names(1:ncol(Z))
sigma.hat = match.arg(sigma.hat)
o <- order(groups)
Z <- Z[,o]
groups <- factor(groups[o])
uniqueGroups <- levels(groups)
ZZ <- lapply(uniqueGroups, function(i) Z[,groups==i, drop=FALSE])
names(ZZ) <- uniqueGroups
sumZ <- sapply(which.mu, function(i) rowSums(ZZ[[i]]))
nGroups = length(uniqueGroups)
sigma2 <- sigma0ld <- rep(ifelse(sigma0>0, sigma0,1), nGroups)
iter = 1
mu <- mu0ld <- rep(0, nGroups)
names(mu) <- uniqueGroups
beta = rep(1,ncol(Z)+length(which.mu))
beta0ld = rep(0,ncol(Z)+length(which.mu))
sigma2.mu = sigma0
nu=0
penalize.mu = FALSE
if(!is.null(which.mu))
if(!penalize.mu)
sumTerm <- "sumZ"
else
sumTerm <- "ridge(sumZ, theta=1/sigma2.mu, scale=FALSE)"
else sumTerm <- character(0)
while((max(abs(beta-beta0ld)) > tol | max(abs(mu - mu0ld)) > tol | max(abs(sigma2 - sigma0ld)) > tol) & iter < max.iter){
beta0ld = beta
sigma0ld <- sigma2
mu0ld <- mu
formula <- formula(paste("surv ~", paste(c(sapply(1:nGroups, function(i) paste("ridge(ZZ[[",i,"]], theta=1/sigma2[",i,"], scale=FALSE)", sep="")),
sumTerm,"strata(strata)"),
collapse=" + ")))
fit <- coxph(formula, ...)
if(any(is.na(coef(fit)))){
warning(paste("NA during estimation (iter: ", iter, ", coef: ", paste(which(is.na(coef(fit)[order(o)])), sep=","), ")", sep=""))
break
}
if(!is.null(which.mu))
mu[which.mu] <- coef(fit)[-(1:ncol(Z))]
if(verbose) cat("mu", mu, "\n", sep="\t")
names(fit$df) <- c(uniqueGroups, rep("Offset", length(which.mu)>0))
if(verbose) cat("df", fit$df,"\n", sep="\t")
sigma2 = sapply(uniqueGroups, function(i){
index <- which(groups==i) #& fit$coefficients > beta.thresh
if(sigma.hat=="BLUP")
(nu * sigma0 + sum((fit$coefficients[index])^2 ))/(nu + length(index))
else if(sigma.hat=="df")
(nu * sigma0 + sum((fit$coefficients[index])^2 ))/(nu + fit$df[i])
else if(sigma.hat == "MLE")
(nu * sigma0 + sum((fit$coefficients[index])^2 ) + sum(diag(solve(solve(fit$var)[index,index]))))/(nu + length(index))
else if(sigma.hat == "REML")
(nu * sigma0 + sum((fit$coefficients[index])^2 ) + sum(diag(fit$var)[index]))/(nu + length(index))
})
if(verbose) {
cat("sigma2", sigma2, "\n", sep="\t")
cat("loglik:", fit$loglik - c(0,fit$penalty[2] + 1/2 * sum(log(sigma2[groups]))),"\n", sep="\t")
}
if(penalize.mu){
if(sigma.hat=="BLUP")
sigma2.mu = (sigma0 * nu + sum((mu-0)^2)) / (nu + length(mu))
else if(sigma.hat=="df")
sigma2.mu = (sigma0 * nu + sum((mu-0)^2)) / (nu + fit$df["Offset"])
else if(sigma.hat == "MLE")
sigma2.mu = (nu * sigma0 + sum((mu-0)^2 ) + sum(diag(solve(solve(fit$var)[-(1:ncol(Z)),-(1:ncol(Z))]))))/(nu + length(mu))
else if(sigma.hat == "REML")
sigma2.mu = (nu * sigma0 + sum((mu-0)^2 ) + sum(diag(fit$var)[-(1:ncol(Z))]))/(nu + length(mu))
}
beta = fit$coefficients
iter = iter+1
}
if(iter == max.iter)
warning("Did not converge after ", max.iter, " iterations.")
fit$iter[1] <- iter
fit$sigma2 = sigma0ld
names(fit$sigma2) <- uniqueGroups
#fit$sigma2.mu = sigma2.mu
fit$mu = mu
fit$Z = data.frame(Z[,order(o)],trans=trans,strata=strata)
fit$surv = surv
fit$tmat=tmat
C <- rbind(diag(1, ncol(Z)),t(as.matrix(MakeInteger(groups)[which.mu]))) ## map from centred to uncentred coefficients
fit$groups = groups[order(o)]
var = fit$var
var2 = fit$var2
colnames(var) <- rownames(var) <- colnames(var2) <- rownames(var2) <- rownames(C) <- c(colnames(Z), which.mu)
colnames(C) <- colnames(Z)
p <- ncol(Z)
i <- c(order(o), (1:ncol(var))[-p:-1])
j <- order(o)
fit$C <- C[i,j]
fit$Hinv <- var[i,i] ## Hinv
fit$V <- var2[i,i] ## Hinv I Hinv
fit$z <- (fit$coefficients / sqrt(diag(var)))[i] ## z-scores of centred coefficients
fit$z2 <- (fit$coefficients / sqrt(diag(var2)))[i] ## z-scores of centred coefficients (var2)
fit$var = (t(C) %*% var %*% C)[j,j] ## covariance of uncentred coef
fit$var2 = (t(C) %*% var2 %*% C)[j,j] ## covariance of uncentred coef (var2)
fit$mu.var = var[-(1:p),-(1:p)] ## covariance of mean
fit$mu.var2 = var2[-(1:p),-(1:p)] ## covariance of mean (var2)
fit$means = fit$means[1:p][j]
fit$coefficients <- (fit$coefficients %*% C)[j]
names(fit$means) <- names(fit$coefficients) <- namesZ
fit$terms <- fit$terms[1:length(uniqueGroups)]
fit$penalized.loglik <- fit$loglik[2] - fit$penalty[2] - 1/2 * sum(log(fit$sigma2[groups]))
## Fake call for predict.coxph and survfit.coxph
call <- match.call()
if(Z.df){
call["data"] <- call["Z"]
formula <- as.formula(paste(as.character(call["surv"]),"~",paste(colnames(Z)[j], collapse="+"),"+strata(strata)"))
}else{
formula <- as.formula(paste(as.character(call["surv"]),"~",as.character(call["Z"])))
}
attr(formula,".Environment") <- parent.frame()
fit$formula <- formula
call["formula"] <- call("foo",formula=formula)["formula"]
fit$terms <- terms(formula)
attr(fit$terms,"specials")<-list(strata=NULL,cluster=NULL)
attr(fit$terms,"specials")$strata<-length(attr(fit$terms,"term.labels"))+1
fit$call <- call
class(fit) <- c("coxrfx", class(fit))
return(fit)
}
#' A summary method for CoxRFX models
#'
#' This function prints the point estimates of parameters and
#' hyperparameters contained in a \code{coxrfx} object.
#'
#'
#' @param object A \code{coxrfx} object
#' (obtained by running the function \code{CoxRFX}).
#' @param ... Further arguments passed to or from other methods.
#' @return NULL
#'
#' @author Rui Costa
#' @export
#' @method summary coxrfx
summary.coxrfx <- function(object,...){
cat("Hyperparameters:\n")
show(format(data.frame(mean=object$mu, variance=object$sigma2), digits=2))
cat("\n Parameters: \n")
show(format(data.frame(group=object$groups, estimate=object$coefficients,row.names = names(object$coefficients)), digits=2))
}
#' Print method for CoxRFX objects
#'
#' This function implicitly calls summary.coxrfx().
#' @param x A \code{coxrfx} object
#' @param ... further arguments passed to or from other methods.
#' @return NULL
#'
#' @author Moritz Gerstung & Rui Costa
#' @export
#' @method print coxrfx
print.coxrfx <- function(x,...){
summary.coxrfx(x)
}
#' Print method for \code{msfit} objects
#' generated by \code{msfit_generic}
#'
#' This method is a simple call to \code{print.default}.
#' Its main purpose is to override \code{print.coxrfx}
#' when printing an
#' object of double class \code{msfit} and \code{coxrfx}.
#' @param x An object of class \code{msfit} or double class \code{msfit}
#' and \code{coxrfx}.
#' @param ... Further arguments passed to or from other methods.
#' @return NULL
#'
#' @author Rui Costa
#' @export
#' @method print msfit
print.msfit <- function(x,...){
#mstate:::summary.msfit(x)
print.default(x,...)
}
#' Convert factor to integer.
#' @param v A factor vector.
#' @return A data.frame with columns corresponding to levels in the factor.
#' @details An internal function of \code{CoxRFX}, not meant
#' to called directly by the user.
#' @author Moritz Gerstung
#' @seealso \code{\link{CoxRFX}}
MakeInteger <- function(v){
res <- as.data.frame(lapply(levels(v), `==`, v))
colnames(res) <- levels(v)
res + 0
}
#coxph function from survival package version 2.44-1.1
#followed by non-exported functions of the same package (potentially used by coxph)
coxph<-function (formula, data, weights, subset, na.action, init, control,
ties = c("efron", "breslow", "exact"), singular.ok = TRUE,
robust = FALSE, model = FALSE, x = FALSE, y = TRUE, tt, method = ties,
...)
{
ties <- match.arg(ties)
Call <- match.call()
extraArgs <- list(...)
if (length(extraArgs)) {
controlargs <- names(formals(coxph.control))
indx <- pmatch(names(extraArgs), controlargs, nomatch = 0L)
if (any(indx == 0L))
stop(gettextf("Argument %s not matched", names(extraArgs)[indx ==
0L]), domain = NA)
}
if (missing(control))
control <- coxph.control(...)
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
names(Call), nomatch = 0)
if (indx[1] == 0)
stop("A formula argument is required")
temp <- Call[c(1, indx)]
temp[[1L]] <- quote(stats::model.frame)
special <- c("strata", "cluster", "tt")
temp$formula <- if (missing(data))
terms(formula, special)
else terms(formula, special, data = data)
if (!is.null(attr(temp$formula, "specials")$tt)) {
coxenv <- new.env(parent = environment(formula))
assign("tt", function(x) x, env = coxenv)
environment(temp$formula) <- coxenv
}
mf <- eval(temp, parent.frame())
if (nrow(mf) == 0)
stop("No (non-missing) observations")
Terms <- terms(mf)
Y <- model.extract(mf, "response")
if (!inherits(Y, "Surv"))
stop("Response must be a survival object")
type <- attr(Y, "type")
if (type != "right" && type != "counting")
stop(paste("Cox model doesn't support \"", type, "\" survival data",
sep = ""))
data.n <- nrow(Y)
if (control$timefix)
Y <- aeqSurv(Y)
if (length(attr(Terms, "variables")) > 2) {
ytemp <- terms.inner(formula[1:2])
xtemp <- terms.inner(formula[-2])
if (any(!is.na(match(xtemp, ytemp))))
warning("a variable appears on both the left and right sides of the formula")
}
strats <- attr(Terms, "specials")$strata
if (length(strats)) {
stemp <- untangle.specials(Terms, "strata", 1)
if (length(stemp$vars) == 1)
strata.keep <- mf[[stemp$vars]]
else strata.keep <- strata(mf[, stemp$vars], shortlabel = TRUE)
strats <- as.numeric(strata.keep)
}
timetrans <- attr(Terms, "specials")$tt
if (missing(tt))
tt <- NULL
if (length(timetrans)) {
timetrans <- untangle.specials(Terms, "tt")
ntrans <- length(timetrans$terms)
if (is.null(tt)) {
tt <- function(x, time, riskset, weights) {
obrien <- function(x) {
r <- rank(x)
(r - 0.5)/(0.5 + length(r) - r)
}
unlist(tapply(x, riskset, obrien))
}
}
if (is.function(tt))
tt <- list(tt)
if (is.list(tt)) {
if (any(!sapply(tt, is.function)))
stop("The tt argument must contain function or list of functions")
if (length(tt) != ntrans) {
if (length(tt) == 1) {
temp <- vector("list", ntrans)
for (i in 1:ntrans) temp[[i]] <- tt[[1]]
tt <- temp
}
else stop("Wrong length for tt argument")
}
}
else stop("The tt argument must contain a function or list of functions")
if (ncol(Y) == 2) {
if (length(strats) == 0) {
sorted <- order(-Y[, 1], Y[, 2])
newstrat <- rep.int(0L, nrow(Y))
newstrat[1] <- 1L
}
else {
sorted <- order(strats, -Y[, 1], Y[, 2])
newstrat <- as.integer(c(1, 1 * (diff(strats[sorted]) !=
0)))
}
if (storage.mode(Y) != "double")
storage.mode(Y) <- "double"
counts <- .Call(Ccoxcount1, Y[sorted, ], as.integer(newstrat))
tindex <- sorted[counts$index]
}
else {
if (length(strats) == 0) {
sort.end <- order(-Y[, 2], Y[, 3])
sort.start <- order(-Y[, 1])
newstrat <- c(1L, rep(0, nrow(Y) - 1))
}
else {
sort.end <- order(strats, -Y[, 2], Y[, 3])
sort.start <- order(strats, -Y[, 1])
newstrat <- c(1L, as.integer(diff(strats[sort.end]) !=
0))
}
if (storage.mode(Y) != "double")
storage.mode(Y) <- "double"
counts <- .Call(Ccoxcount2, Y, as.integer(sort.start -
1L), as.integer(sort.end - 1L), as.integer(newstrat))
tindex <- counts$index
}
Y <- Surv(rep(counts$time, counts$nrisk), counts$status)
type <- "right"
mf <- mf[tindex, ]
strats <- rep(1:length(counts$nrisk), counts$nrisk)
weights <- model.weights(mf)
if (!is.null(weights) && any(!is.finite(weights)))
stop("weights must be finite")
tcall <- attr(Terms, "variables")[timetrans$terms + 2]
pvars <- attr(Terms, "predvars")
pmethod <- sub("makepredictcall.", "", as.vector(methods("makepredictcall")))
for (i in 1:ntrans) {
newtt <- (tt[[i]])(mf[[timetrans$var[i]]], Y[, 1],
strats, weights)
mf[[timetrans$var[i]]] <- newtt
nclass <- class(newtt)
if (any(nclass %in% pmethod)) {
dummy <- as.call(list(as.name(class(newtt)[1]),
tcall[[i]][[2]]))
ptemp <- makepredictcall(newtt, dummy)
pvars[[timetrans$terms[i] + 2]] <- ptemp
}
}
attr(Terms, "predvars") <- pvars
}
cluster <- attr(Terms, "specials")$cluster
if (length(cluster)) {
robust <- TRUE
tempc <- untangle.specials(Terms, "cluster", 1:10)
ord <- attr(Terms, "order")[tempc$terms]
if (any(ord > 1))
stop("Cluster can not be used in an interaction")
cluster <- strata(mf[, tempc$vars], shortlabel = TRUE)
dropterms <- tempc$terms
xlevels <- .getXlevels(Terms[-tempc$terms], mf)
}
else {
if (missing(robust))
robust <- FALSE
xlevels <- .getXlevels(Terms, mf)
dropterms <- NULL
}
contrast.arg <- NULL
attr(Terms, "intercept") <- 1
stemp <- untangle.specials(Terms, "strata", 1)
hasinteractions <- FALSE
if (length(stemp$vars) > 0) {
for (i in stemp$vars) {
if (any(attr(Terms, "order")[attr(Terms, "factors")[i,
] > 0] > 1))
hasinteractions <- TRUE
}
if (!hasinteractions)
dropterms <- c(dropterms, stemp$terms)
}
if (length(dropterms)) {
Terms2 <- Terms[-dropterms]
X <- model.matrix(Terms2, mf, constrasts = contrast.arg)
temp <- attr(X, "assign")
shift <- sort(dropterms)
temp <- temp + 1 * (shift[1] <= temp)
if (length(shift) == 2)
temp + 1 * (shift[2] <= temp)
attr(X, "assign") <- temp
}
else X <- model.matrix(Terms, mf, contrasts = contrast.arg)
if (!all(is.finite(X)))
stop("data contains an infinite predictor")
Xatt <- attributes(X)
if (hasinteractions)
adrop <- c(0, untangle.specials(Terms, "strata")$terms)
else adrop <- 0
xdrop <- Xatt$assign %in% adrop
X <- X[, !xdrop, drop = FALSE]
attr(X, "assign") <- Xatt$assign[!xdrop]
attr(X, "contrasts") <- Xatt$contrasts
offset <- model.offset(mf)
if (is.null(offset) | all(offset == 0))
offset <- rep(0, nrow(mf))
else if (any(!is.finite(offset) | !is.finite(exp(offset))))
stop("offsets must lead to a finite risk score")
weights <- model.weights(mf)
if (!is.null(weights) && any(!is.finite(weights)))
stop("weights must be finite")
assign <- attrassign(X, Terms)
contr.save <- attr(X, "contrasts")
if (missing(init))
init <- NULL
else {
if (length(init) != ncol(X))
stop("wrong length for init argument")
temp <- X %*% init - sum(colMeans(X) * init) + offset
if (any(exp(temp) > .Machine$double.xmax) || all(exp(temp) ==
0))
stop("initial values lead to overflow or underflow of the exp function")
}
if (sum(Y[, ncol(Y)]) == 0) {
ncoef <- ncol(X)
ctemp <- rep(NA, ncoef)
names(ctemp) <- colnames(X)
concordance = c(concordant = 0, discordant = 0, tied.x = 0,
tied.y = 0, tied.xy = 0, concordance = NA, std = NA,
timefix = FALSE)
rval <- list(coefficients = ctemp, var = matrix(0, ncoef,
ncoef), loglik = c(0, 0), score = 0, iter = 0, linear.predictors = offset,
residuals = rep(0, data.n), means = colMeans(X),
method = method, n = data.n, nevent = 0, terms = Terms,
assign = assign, concordance = concordance, y = Y,
call = Call)
class(rval) <- "coxph"
return(rval)
}
pterms <- sapply(mf, inherits, "coxph.penalty")
if (any(pterms)) {
pattr <- lapply(mf[pterms], attributes)
pname <- names(pterms)[pterms]
ord <- attr(Terms, "order")[match(pname, attr(Terms,
"term.labels"))]
if (any(ord > 1))
stop("Penalty terms cannot be in an interaction")
pcols <- assign[match(pname, names(assign))]
fit <- coxpenal.fit(X, Y, strats, offset, init = init,
control, weights = weights, method = method, row.names(mf),
pcols, pattr, assign)
}
else {
if (method == "breslow" || method == "efron") {
if (type == "right")
fitter <- get("coxph.fit")
else fitter <- get("agreg.fit")
}
else if (method == "exact") {
if (type == "right")
fitter <- get("coxexact.fit")
else fitter <- get("agexact.fit")
}
else stop(paste("Unknown method", method))
fit <- fitter(X, Y, strats, offset, init, control, weights = weights,
method = method, row.names(mf))
}
if (is.character(fit)) {
fit <- list(fail = fit)
class(fit) <- "coxph"
}
else {
if (!is.null(fit$coefficients) && any(is.na(fit$coefficients))) {
vars <- (1:length(fit$coefficients))[is.na(fit$coefficients)]
msg <- paste("X matrix deemed to be singular; variable",
paste(vars, collapse = " "))
if (!singular.ok)
stop(msg)
}
fit$n <- data.n
fit$nevent <- sum(Y[, ncol(Y)])
fit$terms <- Terms
fit$assign <- assign
class(fit) <- fit$class
fit$class <- NULL
if (robust && !is.null(fit$coefficients) && !all(is.na(fit$coefficients))) {
fit$naive.var <- fit$var
fit2 <- c(fit, list(x = X, y = Y, weights = weights))
if (length(strats))
fit2$strata <- strats
if (length(cluster)) {
temp <- residuals.coxph(fit2, type = "dfbeta",
collapse = cluster, weighted = TRUE)
if (is.null(init))
fit2$linear.predictors <- 0 * fit$linear.predictors
else fit2$linear.predictors <- c(X %*% init)
temp0 <- residuals.coxph(fit2, type = "score",
collapse = cluster, weighted = TRUE)
}
else {
temp <- residuals.coxph(fit2, type = "dfbeta",
weighted = TRUE)
fit2$linear.predictors <- 0 * fit$linear.predictors
temp0 <- residuals.coxph(fit2, type = "score",
weighted = TRUE)
}
fit$var <- t(temp) %*% temp
u <- apply(as.matrix(temp0), 2, sum)
fit$rscore <- coxph.wtest(t(temp0) %*% temp0, u,
control$toler.chol)$test
}
if (length(fit$coefficients) && is.null(fit$wald.test)) {
nabeta <- !is.na(fit$coefficients)
if (is.null(init))
temp <- fit$coefficients[nabeta]
else temp <- (fit$coefficients - init[1:length(fit$coefficients)])[nabeta]
fit$wald.test <- coxph.wtest(fit$var[nabeta, nabeta],
temp, control$toler.chol)$test
}
if (length(cluster))
temp <- concordancefit(Y, fit$linear.predictors,
strats, weights, cluster = cluster, reverse = TRUE,
timefix = FALSE)
else temp <- concordancefit(Y, fit$linear.predictors,
strats, weights, reverse = TRUE, timefix = FALSE)
if (is.matrix(temp$count))
fit$concordance <- c(colSums(temp$count), concordance = temp$concordance,
std = sqrt(temp$var))
else fit$concordance <- c(temp$count, concordance = temp$concordance,
std = sqrt(temp$var))
na.action <- attr(mf, "na.action")
if (length(na.action))
fit$na.action <- na.action
if (model) {
if (length(timetrans)) {
mf[[".surv."]] <- Y
mf[[".strata."]] <- strats
stop("'model=TRUE' not supported for models with tt terms")
}
fit$model <- mf
}
if (x) {
fit$x <- X
if (length(strats)) {
if (length(timetrans))
fit$strata <- strats
else fit$strata <- strata.keep
}
}
if (y)
fit$y <- Y
}
if (!is.null(weights) && any(weights != 1))
fit$weights <- weights
names(fit$means) <- names(fit$coefficients)
fit$formula <- formula(Terms)
if (length(xlevels) > 0)
fit$xlevels <- xlevels
fit$contrasts <- contr.save
if (any(offset != 0))
fit$offset <- offset
fit$call <- Call
fit
}
terms.inner<-function (x)
{
if (inherits(x, "formula")) {
if (length(x) == 3)
c(terms.inner(x[[2]]), terms.inner(x[[3]]))
else terms.inner(x[[2]])
}
else if (inherits(x, "call") && (x[[1]] != as.name("$") &&
x[[1]] != as.name("["))) {
if (x[[1]] == "+" || x[[1]] == "*" || x[[1]] == "-") {
if (length(x) == 3)
c(terms.inner(x[[2]]), terms.inner(x[[3]]))
else terms.inner(x[[2]])
}
else if (x[[1]] == as.name("Surv"))
unlist(lapply(x[-1], terms.inner))
else terms.inner(x[[2]])
}
else (deparse(x))
}
coxpenal.df<-function (hmat, hinv, fdiag, assign.list, ptype, nvar, pen1,
pen2, sparse)
{
if (ptype == 1 & nvar == 0) {
hdiag <- 1/fdiag
list(fvar2 = (hdiag - pen1) * fdiag^2, df = sum((hdiag -
pen1) * fdiag), fvar = fdiag, trH = sum(fdiag))
}
else if (ptype == 2) {
hmat.full <- t(hmat) %*% (ifelse(fdiag == 0, 0, 1/fdiag) *
hmat)
hinv.full <- hinv %*% (fdiag * t(hinv))
if (length(pen2) == length(hmat.full))
imat <- hmat.full - pen2
else imat <- hmat.full - diag(pen2)
var <- hinv.full %*% imat %*% hinv.full
if (length(assign.list) == 1)
list(var2 = var, df = sum(imat * hinv.full), trH = sum(diag(hinv.full)),
var = hinv.full)
else {
df <- trH <- NULL
d2 <- diag(hinv.full)
for (i in assign.list) {
temp <- coxph.wtest(hinv.full[i, i], var[i, i])$solve
if (is.matrix(temp))
df <- c(df, sum(diag(temp)))
else df <- c(df, sum(temp))
trH <- c(trH, sum(d2[i]))
}
list(var2 = var, df = df, trH = trH, var = hinv.full)
}
}
else {
nf <- length(fdiag) - nvar
nr1 <- 1:nf
nr2 <- (nf + 1):(nf + nvar)
d1 <- fdiag[nr1]
d2 <- fdiag[nr2]
temp <- t(hinv[nr1, ])
temp2 <- t(hinv[nr2, , drop = FALSE])
A.diag <- d1 + c(rep(1, nvar) %*% (temp^2 * d2))
B <- hinv[nr1, ] %*% (d2 * temp2)
C <- hinv[nr2, ] %*% (d2 * temp2)
var2 <- C - t(B) %*% (pen1 * B)
if (ptype == 3) {
hmat.22 <- t(hmat) %*% (ifelse(fdiag == 0, 0, 1/fdiag) *
hmat)
temp <- C - coxph.wtest(hmat.22, diag(nvar))$solve
if (nvar == 1) {
var2 <- var2 - C * pen2 * C
temp2 <- c(temp * pen2)
}
else if (length(pen2) == nvar) {
var2 <- var2 - C %*% (pen2 * C)
temp2 <- sum(diag(temp) * pen2)
}
else {
var2 <- var2 - C %*% matrix(pen2, nvar) %*% C
temp2 <- sum(diag(temp * pen2))
}
}
else temp2 <- 0
df <- trH <- NULL
cdiag <- diag(C)
for (i in 1:length(assign.list)) {
if (sparse == i) {
df <- c(df, nf - (sum(A.diag * pen1) + temp2))
trH <- c(trH, sum(A.diag))
}
else {
j <- assign.list[[i]]
temp <- coxph.wtest(C[j, j], var2[j, j])$solve
if (is.matrix(temp))
df <- c(df, sum(diag(temp)))
else df <- c(df, sum(temp))
trH <- c(trH, sum(cdiag[j]))
}
}
list(var = C, df = df, trH = trH, fvar = A.diag, var2 = var2)
}
}
coxpenal.fit<-function (x, y, strata, offset, init, control, weights, method,
rownames, pcols, pattr, assign)
{
eps <- control$eps
n <- nrow(y)
if (is.matrix(x))
nvar <- ncol(x)
else if (length(x) == 0)
stop("Must have an X variable")
else nvar <- 1
if (missing(offset) || is.null(offset))
offset <- rep(0, n)
if (missing(weights) || is.null(weights))
weights <- rep(1, n)
else {
if (any(weights <= 0))
stop("Invalid weights, must be >0")
}
if (ncol(y) == 3) {
if (length(strata) == 0) {
sorted <- cbind(order(-y[, 2], y[, 3]), order(-y[,
1])) - 1L
newstrat <- as.integer(n)
}
else {
sorted <- cbind(order(strata, -y[, 2], y[, 3]), order(strata,
-y[, 1])) - 1L
newstrat <- as.integer(cumsum(table(strata)))
}
status <- y[, 3]
andersen <- TRUE
}
else {
if (length(strata) == 0) {
sorted <- order(-y[, 1], y[, 2]) - 1L
newstrat <- as.integer(n)
}
else {
sorted <- order(strata, -y[, 1], y[, 2]) - 1L
newstrat <- as.integer(cumsum(table(strata)))
}
status <- y[, 2]
andersen <- FALSE
}
n.eff <- sum(y[, ncol(y)])
npenal <- length(pattr)
if (npenal == 0 || length(pcols) != npenal)
stop("Invalid pcols or pattr arg")
sparse <- sapply(pattr, function(x) !is.null(x$sparse) &&
x$sparse)
if (sum(sparse) > 1)
stop("Only one sparse penalty term allowed")
pterms <- rep(0, length(assign))
names(pterms) <- names(assign)
pindex <- rep(0, npenal)
for (i in 1:npenal) {
temp <- unlist(lapply(assign, function(x, y) (length(x) ==
length(y) && all(x == y)), pcols[[i]]))
if (sparse[i])
pterms[temp] <- 2
else pterms[temp] <- 1
pindex[i] <- (seq(along.with = temp))[temp]
}
if ((sum(pterms == 2) != sum(sparse)) || (sum(pterms > 0) !=
npenal))
stop("pcols and assign arguments disagree")
if (any(pindex != sort(pindex))) {
temp <- order(pindex)
pindex <- pindex[temp]
pcols <- pcols[temp]
pattr <- pattr[temp]
}
ptype <- any(sparse) + 2 * (any(!sparse))
f.expr1 <- function(coef) NULL
f.expr2 <- function(coef) NULL
if (any(sparse)) {
sparse.attr <- (pattr[sparse])[[1]]
fcol <- unlist(pcols[sparse])
if (length(fcol) > 1)
stop("Sparse term must be single column")
xx <- x[, -fcol, drop = FALSE]
for (i in 1:length(assign)) {
j <- assign[[i]]
if (j[1] > fcol)
assign[[i]] <- j - 1
}
for (i in 1:npenal) {
j <- pcols[[i]]
if (j[1] > fcol)
pcols[[i]] <- j - 1
}
frailx <- x[, fcol]
frailx <- match(frailx, sort(unique(frailx)))
nfrail <- max(frailx)
nvar <- nvar - 1
pfun1 <- sparse.attr$pfun
f.expr1 <- function(coef) {
coxlist1$coef <- coef
if (is.null(extra1))
temp <- pfun1(coef, theta1, n.eff)
else temp <- pfun1(coef, theta1, n.eff, extra1)
if (!is.null(temp$recenter))
coxlist1$coef <- coxlist1$coef - as.double(temp$recenter)
if (!temp$flag) {
coxlist1$first <- -as.double(temp$first)
coxlist1$second <- as.double(temp$second)
}
coxlist1$penalty <- -as.double(temp$penalty)
coxlist1$flag <- as.logical(temp$flag)
if (any(sapply(coxlist1, length) != c(rep(nfrail,
3), 1, 1)))
stop("Incorrect length in coxlist1")
coxlist1
}
if (!is.null(getOption("survdebug")))
debug(f.expr1)
coxlist1 <- list(coef = double(nfrail), first = double(nfrail),
second = double(nfrail), penalty = 0, flag = FALSE)
}
else {
xx <- x
frailx <- 0
nfrail <- 0
}
if (sum(!sparse) > 0) {
full.imat <- !all(unlist(lapply(pattr, function(x) x$diag)))
ipenal <- (1:length(pattr))[!sparse]
f.expr2 <- function(coef) {
coxlist2$coef <- coef
pentot <- 0
for (i in ipenal) {
pen.col <- pcols[[i]]
coef <- coxlist2$coef[pen.col]
if (is.null(extralist[[i]]))
temp <- ((pattr[[i]])$pfun)(coef, thetalist[[i]],
n.eff)
else temp <- ((pattr[[i]])$pfun)(coef, thetalist[[i]],
n.eff, extralist[[i]])
if (!is.null(temp$recenter))
coxlist2$coef[pen.col] <- coxlist2$coef[pen.col] -
temp$recenter
if (temp$flag)
coxlist2$flag[pen.col] <- TRUE
else {
coxlist2$flag[pen.col] <- FALSE
coxlist2$first[pen.col] <- -temp$first
if (full.imat) {
tmat <- matrix(coxlist2$second, nvar, nvar)
tmat[pen.col, pen.col] <- temp$second
coxlist2$second <- c(tmat)
}
else coxlist2$second[pen.col] <- temp$second
}
pentot <- pentot - temp$penalty
}
coxlist2$penalty <- as.double(pentot)
if (any(sapply(coxlist2, length) != length2))
stop("Length error in coxlist2")
coxlist2
}
if (!is.null(getOption("survdebug")))
debug(f.expr2)
if (full.imat) {
coxlist2 <- list(coef = double(nvar), first = double(nvar),
second = double(nvar * nvar), penalty = 0, flag = rep(FALSE,
nvar))
length2 <- c(nvar, nvar, nvar * nvar, 1, nvar)
}
else {
coxlist2 <- list(coef = double(nvar), first = double(nvar),
second = double(nvar), penalty = 0, flag = rep(FALSE,
nvar))
length2 <- c(nvar, nvar, nvar, 1, nvar)
}
}
else full.imat <- FALSE
if (nfrail > 0)
finit <- rep(0, nfrail)
else finit <- 0
if (!missing(init) && !is.null(init)) {
if (length(init) != nvar) {
if (length(init) == (nvar + nfrail)) {
finit <- init[-(1:nvar)]
init <- init[1:nvar]
}
else stop("Wrong length for inital values")
}
}
else init <- double(nvar)
cfun <- lapply(pattr, function(x) x$cfun)
parmlist <- lapply(pattr, function(x, eps) c(x$cparm, eps2 = eps),
sqrt(eps))
extralist <- lapply(pattr, function(x) x$pparm)
iterlist <- vector("list", length(cfun))
thetalist <- vector("list", length(cfun))
printfun <- lapply(pattr, function(x) x$printfun)
for (i in 1:length(cfun)) {
temp <- (cfun[[i]])(parmlist[[i]], iter = 0)
if (sparse[i]) {
theta1 <- temp$theta
extra1 <- extralist[[i]]
}
thetalist[[i]] <- temp$theta
iterlist[[i]] <- temp
}
temp1 <- c("x", "coef", "plik", "loglik", "status", "neff",
"df", "trH")
temp2 <- c("frailx", "coxfit$fcoef", "loglik1", "coxfit$loglik",
"status", "n.eff")
temp3 <- c("xx[,pen.col]", "coxfit$coef[pen.col]", "loglik1",
"coxfit$loglik", "status", "n.eff")
calls <- vector("expression", length(cfun))
cargs <- lapply(pattr, function(x) x$cargs)
for (i in 1:length(cfun)) {
tempchar <- paste("(cfun[[", i, "]])(parmlist[[", i,
"]], iter,", "iterlist[[", i, "]]")
temp2b <- c(temp2, paste("pdf[", i, "]"), paste("trH[",
i, "]"))
temp3b <- c(temp3, paste("pdf[", i, "]"), paste("trH[",
i, "]"))
if (length(cargs[[i]]) == 0)
calls[i] <- parse(text = paste(tempchar, ")"))
else {
temp <- match(cargs[[i]], temp1)
if (any(is.na(temp)))
stop(paste((cargs[[i]])[is.na(temp)], "not matched"))
if (sparse[i])
temp4 <- paste(temp2b[temp], collapse = ",")
else temp4 <- paste(temp3b[temp], collapse = ",")
calls[i] <- parse(text = paste(paste(tempchar, temp4,
sep = ","), ")"))
}
}
need.df <- any(!is.na(match(c("df", "trH"), unlist(cargs))))
varnames <- dimnames(xx)[[2]]
for (i in 1:length(cfun)) {
if (!is.null(pattr[[i]]$varname))
varnames[pcols[[i]]] <- pattr[[i]]$varname
}
rho <- environment()
storage.mode(y) <- storage.mode(weights) <- "double"
storage.mode(xx) <- storage.mode(offset) <- "double"
if (andersen)
coxfit <- .C(Cagfit5a, as.integer(n), as.integer(nvar),
y, xx, offset, weights, newstrat, sorted, means = double(nvar),
coef = as.double(init), u = double(nvar), loglik = double(1),
as.integer(method == "efron"), as.integer(ptype),
as.integer(full.imat), as.integer(nfrail), as.integer(frailx),
f.expr1, f.expr2, rho)
else coxfit <- .C(Ccoxfit5a, as.integer(n), as.integer(nvar),
y, xx, offset, weights, newstrat, sorted, means = double(nvar),
coef = as.double(init), u = double(nvar), loglik = double(1),
as.integer(method == "efron"), as.integer(ptype), as.integer(full.imat),
as.integer(nfrail), as.integer(frailx), f.expr1, f.expr2,
rho)
loglik0 <- coxfit$loglik
means <- coxfit$means
iter2 <- 0
iterfail <- NULL
thetasave <- unlist(thetalist)
for (outer in 1:control$outer.max) {
if (andersen)
coxfit <- .C(Cagfit5b, iter = as.integer(control$iter.max),
as.integer(n), as.integer(nvar), as.integer(newstrat),
coef = as.double(init), u = double(nvar + nfrail),
hmat = double(nvar * (nvar + nfrail)), hinv = double(nvar *
(nvar + nfrail)), loglik = double(1), flag = integer(1),
as.double(control$eps), as.double(control$toler.chol),
as.integer(method == "efron"), as.integer(nfrail),
fcoef = as.double(finit), fdiag = double(nfrail +
nvar), f.expr1, f.expr2, rho)
else coxfit <- .C(Ccoxfit5b, iter = as.integer(control$iter.max),
as.integer(n), as.integer(nvar), as.integer(newstrat),
coef = as.double(init), u = double(nvar + nfrail),
hmat = double(nvar * (nvar + nfrail)), hinv = double(nvar *
(nvar + nfrail)), loglik = double(1), flag = integer(1),
as.double(control$eps), as.double(control$toler.chol),
as.integer(method == "efron"), as.integer(nfrail),
fcoef = as.double(finit), fdiag = double(nfrail +
nvar), f.expr1, f.expr2, rho)
iter <- outer
iter2 <- iter2 + coxfit$iter
if (coxfit$iter >= control$iter.max)
iterfail <- c(iterfail, iter)
temp <- rep(FALSE, nvar + nfrail)
if (nfrail > 0)
temp[1:nfrail] <- coxlist1$flag
if (ptype > 1)
temp[nfrail + 1:nvar] <- coxlist2$flag
fdiag <- ifelse(temp, 0, coxfit$fdiag)
if (need.df) {
if (nfrail > 0)
temp1 <- coxlist1$second
else temp1 <- 0
if (ptype > 1)
temp2 <- coxlist2$second
else temp2 <- 0
dftemp <- coxpenal.df(matrix(coxfit$hmat, ncol = nvar),
matrix(coxfit$hinv, ncol = nvar), fdiag, assign,
ptype, nvar, temp1, temp2, pindex[sparse])
df <- dftemp$df
var <- dftemp$var
var2 <- dftemp$var2
pdf <- df[pterms > 0]
trH <- dftemp$trH[pterms > 0]
}
if (nfrail > 0)
penalty <- -coxlist1$penalty
else penalty <- 0
if (ptype > 1)
penalty <- penalty - coxlist2$penalty
loglik1 <- coxfit$loglik + penalty
if (iter == 1)
penalty0 <- penalty
done <- TRUE
for (i in 1:length(cfun)) {
pen.col <- pcols[[i]]
temp <- eval(calls[i])
if (sparse[i])
theta1 <- temp$theta
thetalist[[i]] <- temp$theta
iterlist[[i]] <- temp
done <- done & temp$done
}
if (done)
break
if (iter == 1) {
init <- coefsave <- coxfit$coef
finit <- fsave <- coxfit$fcoef
thetasave <- cbind(thetasave, unlist(thetalist))
}
else {
temp <- as.vector(unlist(thetalist))
coefsave <- cbind(coefsave, coxfit$coef)
fsave <- cbind(fsave, coxfit$fcoef)
howclose <- apply((thetasave - temp)^2, 2, sum)
which <- min((1:iter)[howclose == min(howclose)])
if (nvar > 0)
init <- coefsave[, which]
if (nfrail > 0)
finit <- fsave[, which]
thetasave <- cbind(thetasave, temp)
}
}
if (nfrail > 0) {
lp <- offset + coxfit$fcoef[frailx]
if (nvar > 0)
lp <- lp + x[, -fcol, drop = FALSE] %*% coxfit$coef -
sum(means * coxfit$coef)
}
else lp <- offset + as.vector(x %*% coxfit$coef) - sum(means *
coxfit$coef)
if (andersen) {
.C(Cagfit5c, as.integer(nvar))
resid <- .Call(Cagmart3, y, exp(lp), weights, newstrat,
sorted, as.integer(method == "efron"))
}
else {
expect <- .C(Ccoxfit5c, as.integer(n), as.integer(nvar),
as.integer(newstrat), as.integer(method == "efron"),
expect = double(n))$expect
resid <- status - expect
}
names(resid) <- rownames
if (!need.df) {
if (nfrail > 0)
temp1 <- coxlist1$second
else temp1 <- 0
if (ptype > 1)
temp2 <- coxlist2$second
else temp2 <- 0
dftemp <- coxpenal.df(matrix(coxfit$hmat, ncol = nvar),
matrix(coxfit$hinv, ncol = nvar), fdiag, assign,
ptype, nvar, temp1, temp2, pindex[sparse])
df <- dftemp$df
trH <- dftemp$trH
var <- dftemp$var
var2 <- dftemp$var2
}
if (control$iter.max > 1 && length(iterfail) > 0)
warning(paste("Inner loop failed to coverge for iterations",
paste(iterfail, collapse = " ")))
which.sing <- (fdiag[nfrail + 1:nvar] == 0)
coef <- coxfit$coef
names(coef) <- varnames
coef[which.sing] <- NA
names(iterlist) <- names(pterms[pterms > 0])
if (nfrail > 0) {
if (nvar > 0) {
list(coefficients = coef, var = var, var2 = var2,
loglik = c(loglik0, loglik1), iter = c(iter,
iter2), linear.predictors = as.vector(lp),
residuals = resid, means = means, method = method,
class = c("coxph.penal", "coxph"), frail = coxfit$fcoef,
fvar = dftemp$fvar, df = df, df2 = dftemp$df2,
penalty = c(penalty0, penalty), pterms = pterms,
assign2 = assign, history = iterlist, coxlist1 = coxlist1,
printfun = printfun)
}
else {
list(loglik = c(loglik0, loglik1), iter = c(iter,
iter2), linear.predictors = as.vector(lp), residuals = resid,
means = means, method = method, class = c("coxph.penal",
"coxph"), frail = coxfit$fcoef, fvar = dftemp$fvar,
df = df, df2 = dftemp$df2, penalty = c(penalty0,
penalty), pterms = pterms, assign2 = assign,
history = iterlist, printfun = printfun)
}
}
else {
list(coefficients = coef, var = var, var2 = var2, loglik = c(loglik0,
loglik1), iter = c(iter, iter2), linear.predictors = lp,
residuals = resid, means = means, method = method,
class = c("coxph.penal", "coxph"), df = df, df2 = dftemp$df2,
penalty = c(penalty0, penalty), pterms = pterms,
assign2 = assign, history = iterlist, coxlist2 = coxlist2,
printfun = printfun)
}
}
parsecovar1<-function (flist, statedatanames)
{
if (any(sapply(flist, function(x) !inherits(x, "formula"))))
stop("flist must be a list of formulas")
if (any(sapply(flist, length) != 3))
stop("all formulas must have a left and right side")
lhs <- lapply(flist, function(x) x[-3])
rhs <- lapply(flist, function(x) x[-2])
rh2 <- lapply(rhs, function(form) {
parts <- strsplit(deparse(form, width.cutoff = 300, control = NULL),
"/", fixed = TRUE)[[1]]
if (length(parts) == 1) {
ival <- NULL
common <- FALSE
fixed <- FALSE
clear <- FALSE
}
else {
optterms <- terms(formula(paste("~", parts[2])))
ff <- rownames(attr(optterms, "factors"))
index <- match(ff, c("common", "fixed", "init", "clear"))
if (any(is.na(index)))
stop("option not recognized in a covariates formula: ",
paste(ff[is.na(index)], collapse = ", "))
common <- any(index == 1)
fixed <- any(index == 2)
clear <- any(index == 3)
if (any(index == 3)) {
optatt <- attributes(optterms)
j <- optatt$variables[1 + which(index == 3)]
j[[1]] <- as.name("list")
ival <- unlist(eval(j, parent.frame()))
}
else ival <- NULL
}
form <- formula(paste("~ -1 +", parts[1]))
list(common = common, fixed = fixed, clear = clear, ival = ival,
formula = form)
})
pcut <- function(form) {
if (length(form) == 3) {
if (form[[1]] == "+")
c(pcut(form[[2]]), pcut(form[[3]]))
else if (form[[1]] == "~")
pcut(form[[2]])
else list(form)
}
else list(form)
}
lcut <- lapply(lhs, function(x) pcut(x[[2]]))
env1 <- new.env(parent = parent.frame(2))
env2 <- new.env(parent = env1)
if (missing(statedatanames)) {
assign("state", function(...) list(stateid = "state",
values = c(...)), env1)
assign("state", list(stateid = "state"))
}
else {
for (i in statedatanames) {
assign(i, eval(list(stateid = i)), env2)
tfun <- eval(parse(text = paste0("function(...) list(stateid='",
i, "', values=c(...)")))
assign(i, tfun, env1)
}
}
lterm <- lapply(lcut, function(x) {
lapply(x, function(z) {
if (length(z) == 1) {
temp <- eval(z, envir = env2)
if (is.list(temp) && names(temp)[[1]] == "stateid")
temp
else temp
}
else if (length(z) == 3 && z[[1]] == ":")
list(left = eval(z[[2]], envir = env2), right = eval(z[[3]],
envir = env2))
else stop("invalid term: ", deparse(z))
})
})
list(rhs = rh2, lhs = lterm)
}
parsecovar2<-function (covar1, statedata, dformula, Terms, transitions, states)
{
if (is.null(statedata))
statedata <- data.frame(state = states)
else {
if (is.null(statedata$state))
stop("the statedata data set must contain a variable 'state'")
indx1 <- match(states, statedata$state, nomatch = 0)
if (any(indx1 == 0))
stop("statedata does not contain all the possible states: ",
states[indx1 == 0])
statedata <- statedata[indx1, ]
}
allterm <- attr(Terms, "term.labels")
nterm <- length(allterm)
nstate <- length(states)
tmap <- array(0L, dim = c(nterm + 1, nstate, nstate))
dterms <- match(attr(terms.formula(dformula), "term.labels"),
allterm)
dterms <- c(1L, 1L + dterms)
k <- seq(along = dterms)
for (i in 1:nstate) {
for (j in 1:nstate) {
tmap[dterms, i, j] <- k
k <- k + length(k)
}
}
ncoef <- max(tmap)
inits <- NULL
if (!is.null(covar1)) {
for (i in 1:length(covar1$rhs)) {
rhs <- covar1$rhs[[i]]
lhs <- covar1$lhs[[i]]
rterm <- terms.formula(rhs$formula)
rindex <- 1L + match(attr(rterm, "term.labels"),
allterm, nomatch = 0)
if (any(rindex == 1L))
stop("dterm mismatch bug 2")
if (attr(rterm, "intercept") == 1)
rindex <- c(1L, rindex)
state1 <- state2 <- NULL
for (x in lhs) {
if (is.null(x$left))
stop("term found without a :", x)
if (!is.list(x$left) && length(x$left) == 1 &
x$left == 1)
temp1 <- 1:nrow(statedata)
else if (is.list(x$left) && names(x$left)[1] ==
"stateid") {
if (is.null(x$left$value))
stop("state variable with no list of values: ",
x$left$stateid)
else temp1 <- which(statedata[[x$left$stateid]] %in%
x$left$value)
}
else temp1 <- which(statedata$state %in% x$left)
if (!is.list(x$right) && length(x$right) == 1 &&
x$right == 1)
temp2 <- 1:nrow(statedata)
else if (is.list(x$right) && names(x$right)[1] ==
"stateid") {
if (is.null(x$right$value))
stop("state variable with no list of values: ",
x$right$stateid)
else temp2 <- which(statedata[[x$right$stateid]] %in%
x$right$value)
}
else temp2 <- which(statedata$state %in% x$right)
state1 <- c(state1, rep(temp1, length(temp2)))
state2 <- c(state2, rep(temp2, each = length(temp1)))
}
if (rhs$clear)
tmap[-1, state1, state2] <- 0
if (length(rhs$ival))
inits <- c(inits, list(term = rindex, state1 = state1,
state2 = state2, init = rhs$ival))
if (rhs$common)
j <- ncoef + seq_len(length(rindex))
else j <- ncoef + seq_len(length(rindex) * length(state1))
tmap[rindex, state1, state2] <- j
ncoef <- max(j)
}
}
t2 <- transitions[, -1, drop = FALSE]
indx1 <- match(rownames(t2), states)
indx2 <- match(colnames(t2), states)
tmap2 <- matrix(0L, nrow = 1 + nterm, ncol = sum(t2 > 0))
trow <- row(t2)[t2 > 0]
tcol <- col(t2)[t2 > 0]
for (i in 1:length(trow)) tmap2[, i] <- tmap[, indx1[trow[i]],
indx2[tcol[i]]]
dimnames(tmap2) <- list(c("Intercept", allterm), paste(indx1[trow],
indx2[tcol], sep = ":"))
list(tmap = tmap2, inits = inits, mapid = cbind(indx1[trow],
indx2[tcol]))
}
parsecovar3<-function (tmap, Xcol, Xassign)
{
hasintercept <- (Xassign[1] == 0)
cmap <- matrix(0L, length(Xcol) + !hasintercept, ncol(tmap))
cmap[1, ] <- match(tmap[1, ], sort(c(0, unique(tmap[1, ])))) -
1L
xcount <- table(factor(Xassign, levels = 1:max(Xassign)))
mult <- 1 + max(xcount)
j <- 1
for (i in 2:nrow(tmap)) {
k <- seq_len(xcount[i - 1])
cmap[j + k, ] <- ifelse(tmap[i, ] == 0, 0, tmap[i, ] *
mult + rep(k, ncol(tmap)))
j <- j + max(k)
}
cmap[-1, ] <- match(cmap[-1, ], sort(unique(c(0L, cmap[-1,
])))) - 1L
colnames(cmap) <- colnames(tmap)
if (hasintercept)
rownames(cmap) <- Xcol
else rownames(cmap) <- c("(Intercept)", Xcol)
cmap
}
residuals.coxph<-function (object, type = c("martingale", "deviance", "score",
"schoenfeld", "dfbeta", "dfbetas", "scaledsch", "partial"),
collapse = FALSE, weighted = FALSE, ...)
{
type <- match.arg(type)
otype <- type
if (type == "dfbeta" || type == "dfbetas") {
otype <- type
type <- "score"
if (missing(weighted))
weighted <- TRUE
}
if (type == "scaledsch")
type <- "schoenfeld"
n <- length(object$residuals)
rr <- object$residuals
y <- object$y
x <- object[["x"]]
vv <- drop(object$naive.var)
if (is.null(vv))
vv <- drop(object$var)
weights <- object$weights
if (is.null(weights))
weights <- rep(1, n)
strat <- object$strata
method <- object$method
if (method == "exact" && (type == "score" || type == "schoenfeld"))
stop(paste(otype, "residuals are not available for the exact method"))
if (type == "martingale" || type == "partial")
rr <- object$residuals
else {
Terms <- object$terms
if (!inherits(Terms, "terms"))
stop("invalid terms component of object")
strats <- attr(Terms, "specials")$strata
if (is.null(y) || (is.null(x) && type != "deviance")) {
temp <- coxph.getdata(object, y = TRUE, x = TRUE,
stratax = TRUE)
y <- temp$y
x <- temp$x
if (length(strats))
strat <- temp$strata
}
ny <- ncol(y)
status <- y[, ny, drop = TRUE]
if (type != "deviance") {
nstrat <- as.numeric(strat)
nvar <- ncol(x)
if (is.null(strat)) {
ord <- order(y[, ny - 1], -status)
newstrat <- rep(0, n)
}
else {
ord <- order(nstrat, y[, ny - 1], -status)
newstrat <- c(diff(as.numeric(nstrat[ord])) !=
0, 1)
}
newstrat[n] <- 1
x <- x[ord, ]
y <- y[ord, ]
score <- exp(object$linear.predictors)[ord]
}
}
if (type == "schoenfeld") {
if (ny == 2) {
mintime <- min(y[, 1])
if (mintime < 0)
y <- cbind(2 * mintime - 1, y)
else y <- cbind(-1, y)
}
temp <- .C(Ccoxscho, n = as.integer(n), as.integer(nvar),
as.double(y), resid = as.double(x), as.double(score *
weights[ord]), as.integer(newstrat), as.integer(method ==
"efron"), double(3 * nvar))
deaths <- y[, 3] == 1
if (nvar == 1)
rr <- temp$resid[deaths]
else rr <- matrix(temp$resid[deaths], ncol = nvar)
if (weighted)
rr <- rr * weights[deaths]
if (length(strats))
attr(rr, "strata") <- table((strat[ord])[deaths])
time <- c(y[deaths, 2])
if (is.matrix(rr))
dimnames(rr) <- list(time, names(object$coefficients))
else names(rr) <- time
if (otype == "scaledsch") {
ndead <- sum(deaths)
coef <- ifelse(is.na(object$coefficients), 0, object$coefficients)
rr <- drop(rr %*% vv * ndead + rep(coef, each = nrow(rr)))
}
return(rr)
}
if (type == "score") {
if (ny == 2) {
resid <- .C(Ccoxscore, as.integer(n), as.integer(nvar),
as.double(y), x = as.double(x), as.integer(newstrat),
as.double(score), as.double(weights[ord]), as.integer(method ==
"efron"), resid = double(n * nvar), double(2 *
nvar))$resid
}
else {
resid <- .C(Cagscore, as.integer(n), as.integer(nvar),
as.double(y), as.double(x), as.integer(newstrat),
as.double(score), as.double(weights[ord]), as.integer(method ==
"efron"), resid = double(n * nvar), double(nvar *
6))$resid
}
if (nvar > 1) {
rr <- matrix(0, n, nvar)
rr[ord, ] <- matrix(resid, ncol = nvar)
dimnames(rr) <- list(names(object$residuals), names(object$coefficients))
}
else rr[ord] <- resid
if (otype == "dfbeta") {
if (is.matrix(rr))
rr <- rr %*% vv
else rr <- rr * vv
}
else if (otype == "dfbetas") {
if (is.matrix(rr))
rr <- (rr %*% vv) %*% diag(sqrt(1/diag(vv)))
else rr <- rr * sqrt(vv)
}
}
if (weighted)
rr <- rr * weights
if (!is.null(object$na.action)) {
rr <- naresid(object$na.action, rr)
if (is.matrix(rr))
n <- nrow(rr)
else n <- length(rr)
if (type == "deviance")
status <- naresid(object$na.action, status)
}
if (type == "partial") {
rr <- rr + predict(object, type = "terms")
}
if (!missing(collapse)) {
if (length(collapse) != n)
stop("Wrong length for 'collapse'")
rr <- drop(rowsum(rr, collapse))
if (type == "deviance")
status <- drop(rowsum(status, collapse))
}
if (type == "deviance")
sign(rr) * sqrt(-2 * (rr + ifelse(status == 0, 0, status *
log(status - rr))))
else rr
}
coxph.getdata<-function (fit, y = TRUE, x = TRUE, stratax = TRUE, offset = FALSE)
{
ty <- fit[["y"]]
tx <- fit[["x"]]
strat <- fit$strata
Terms <- fit$terms
if (is.null(attr(Terms, "offset")))
offset <- FALSE
if (offset)
x <- TRUE
if (!inherits(Terms, "terms"))
stop("invalid terms component of fit")
strats <- attr(Terms, "specials")$strata
if (length(strats) == 0)
stratax <- FALSE
if ((y && is.null(ty)) || (x && is.null(tx)) || (stratax &&
is.null(strat)) || offset) {
m <- stats::model.frame(fit)
if (y && is.null(ty))
ty <- model.extract(m, "response")
if (offset)
toff <- model.extract(m, "offset")
if ((x || stratax) && is.null(tx)) {
if (stratax) {
temp <- untangle.specials(Terms, "strata", 1)
strat <- strata(m[temp$vars], shortlabel = TRUE)
}
if (x)
tx <- model.matrix(fit, data = m)
}
}
else if (offset)
toff <- fit$linear.predictors - (c(tx %*% fit$coef) -
sum(fit$means * fit$coef))
temp <- list()
if (y)
temp$y <- ty
if (x)
temp$x <- tx
if (stratax)
temp$strata <- strat
if (offset)
temp$offset <- toff
temp
}
# Create a strata variable, possibly from many objects
#
strata <- function(..., na.group=FALSE, shortlabel, sep=', ') {
# First, grab a copy of the call, which will be used to manufacture
# labels for unlabeled arguments
# Then get the arguments as a list
words <- as.character((match.call())[-1])
allf <- list(...)
# If there is only one argument, and it itself is a list, use
# it instead
if(length(allf) == 1 && is.list(ttt <- unclass(allf[[1]]))) allf <- ttt
nterms <- length(allf)
# Keep the names of named args as their label, what was typed otherwise
if (is.null(names(allf))) {
argname <- words[1:nterms]
if (missing(shortlabel))
shortlabel <- all(sapply(allf,
function(x) is.character(x) | inherits(x, 'factor')))
}
else {
argname <- ifelse(names(allf) == '', words[1:nterms], names(allf))
if (missing(shortlabel)) shortlabel <- FALSE
}
# If the arguments are not all the same length, stop now
# Mostly this is to stop calls with an improper object
arglength <- sapply(allf, length)
if (any(arglength != arglength[1]))
stop("all arguments must be the same length")
if (!all(sapply(allf, is.atomic))) stop("all arguments must be vectors")
# Process the first argument
what <- allf[[1]]
if(is.null(levels(what)))
what <- factor(what)
levs <- unclass(what) - 1
wlab <- levels(what)
if (na.group && any(is.na(what))){
# add "NA" as a level
levs[is.na(levs)] <- length(wlab)
wlab <- c(wlab, "NA")
}
if (shortlabel) labs <- wlab
else labs <- paste(argname[1], wlab, sep='=')
# Now march through the other variables, if any
for (i in (1:nterms)[-1]) {
what <- allf[[i]]
if(is.null(levels(what)))
what <- factor(what)
wlab <- levels(what)
wlev <- unclass(what) - 1
if (na.group && any(is.na(wlev))){
wlev[is.na(wlev)] <- length(wlab)
wlab <- c(wlab, "NA")
}
if (!shortlabel) wlab <- format(paste(argname[i], wlab, sep='='))
levs <- wlev + levs*(length(wlab))
labs <- paste(rep(labs, rep(length(wlab), length(labs))),
rep(wlab, length(labs)), sep=sep)
}
levs <- levs + 1
ulevs <- sort(unique(levs[!is.na(levs)]))
levs <- match(levs, ulevs)
labs <- labs[ulevs]
factor(levs, labels=labs)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.