version 3.8-3

Author: Terry M Therneau

DESCRIPTION

@@ -1,8 +1,8 @@
 Title: Survival Analysis
 Priority: recommended
 Package: survival
-Version: 3.7-0
-Date: 2024-06-01
+Version: 3.8-3
+Date: 2024-12-17
 Depends: R (>= 3.5.0)
 Imports: graphics, Matrix, methods, splines, stats, utils
 LazyData: Yes
@@ -22,12 +22,12 @@ Description: Contains the core survival analysis routines, including
 License: LGPL (>= 2)
 URL: https://github.com/therneau/survival
 NeedsCompilation: yes
-Packaged: 2024-06-03 15:17:04 UTC; therneau
+Packaged: 2024-12-17 16:37:18 UTC; therneau
 Author: Terry M Therneau [aut, cre],
   Thomas Lumley [ctb, trl] (original S->R port and R maintainer until
     2009),
   Atkinson Elizabeth [ctb],
   Crowson Cynthia [ctb]
 Maintainer: Terry M Therneau <therneau.terry@mayo.edu>
 Repository: CRAN
-Date/Publication: 2024-06-05 16:30:02 UTC
+Date/Publication: 2024-12-17 20:20:02 UTC

MD5

@@ -35,7 +35,7 @@ export(Surv, Surv2, Surv2data, aeqSurv, aareg, agreg.fit,
        survobrien,  
        survpenal.fit, 
        survreg, survreg.control, survreg.fit,
-       survreg.distributions, survregDtest, tcut,
+       survreg.distributions, survregDtest, tcut, totimeline, fromtimeline,
        tmerge, untangle.specials, yates, yates_setup,
        survConcordance, survConcordance.fit, survfitcoxph.fit)
 S3method('[', aareg)

NAMESPACE

@@ -18,6 +18,12 @@ aareg <- function(formula, data, weights, subset, na.action,
     #  a local copy, doing otherwise messes up future use of update() on
     #  the model object for a user stuck in "+ cluster()" mode.
     if (missing(formula)) stop("a formula argument is required")
+    # make Surv(), strata() resolve to the survival namespace
+    newform <- removeDoubleColonSurv(formula)
+    if (!is.null(newform)) {
+        formula <- newform$formula
+        if (newform$newcall) Call$formula <- formula
+    }
 
     ss <- c("cluster", "offset")
     Terms <- if (missing(data)) terms(formula, specials=ss) else

R/aareg.R

@@ -1,5 +1,13 @@
-# Automatically generated from the noweb directory
 agsurv <- function(y, x, wt, risk, survtype, vartype) {
+    # This is called by coxsurv.fit, once per stratum, to compute all the
+    # peices of a coxph survival curve
+    # y: Surv object
+    # x: covariate matrix
+    # wt: case weights
+    # risk: exp(X beta) = risk score for each observation
+    # survtype =1 Kalbfleisch-Prentice, 2= exp(cum hazard): Breslow,
+    #  3= exp(cum hazard) : Efron
+    #
     nvar <- ncol(as.matrix(x))
     status <- y[,ncol(y)]
     dtime <- y[,ncol(y) -1]
@@ -74,7 +82,13 @@ agsurv <- function(y, x, wt, risk, survtype, vartype) {
     result <- list(n= nrow(y), time=time, n.event=nevent, n.risk=irisk, 
                    n.censor=ncens, hazard=haz, 
                    cumhaz=cumsum(haz), varhaz=varhaz, ndeath=ndeath, 
-                   xbar=apply(matrix(xbar, ncol=nvar),2, cumsum))
-    if (survtype==1) result$surv <- km$inc
+                   xbar= matrix(xbar, ncol=nvar))
+    # varhaz is the first part of the variance; the increment to var(\Lambda)
+    #  at each time, that one we would have if beta were a fixed constant.
+    # Each row of xbar is (mean of those at risk) * lambda; we
+    # need it for the second part of the variance.  The other elements are
+    # what you would expect.
+
+    if (survtype==1) result$surv <- km$inc  # Kalbfleisch-Prentice
     result
 }

R/agsurv.R

@@ -14,18 +14,25 @@ anova.coxphlist <- function (object, test =  'Chisq' ,...) {
         stop("all models must have the same ties option")
 
     responses <- as.character(unlist(lapply(object, 
-				     function(x) deparse(formula(x)[[2]]))))
+				     function(x) deparse1(formula(x)[[2]]))))
     sameresp <- (responses == responses[1])
     if (!all(sameresp)) {
         object <- object[sameresp]
-        warning(paste("Models with response", deparse(responses[!sameresp]), 
+        warning(paste("Models with response", deparse1(responses[!sameresp]), 
             "removed because response differs from", "model 1"))
     }
 
     ns <- sapply(object, function(x) length(x$residuals))
     if (any(ns != ns[1])) 
         stop("models were not all fit to the same size of dataset")
 
+    # verify that all of them used the same strata, if present
+    stemp <- lapply(object, function(x) 
+        untangle.specials(x$terms, "strata")[["vars"]])
+    fail <- (any(sapply(stemp, length) > 0) &&
+             any(sapply(stemp, function(x) !identical(x, stemp[[1]]))))
+    if (fail) stop("models do not have the same strata")
+
     nmodels <- length(object)
     if (nmodels == 1) # only one model remains
         return(anova.coxph(object[[1]], test = test))

R/anova.coxphlist.R

@@ -0,0 +1,94 @@
+# This is has arguments like anova.coxph: a call with a single model is
+#  not allowed.  So: we avoid having separate anova.coxphms and anova.coxphmslist
+#  functions.
+# 
+anova.coxphms <- function (object, ...,  test = c("score", "Wald", "PL")) {
+    if (!inherits(object, "coxphms"))
+        stop ("argument must be the fit of a multistate hazard  model")
+
+    # All the ... args need to be coxphms fits. If any of them
+    #  have a name attached, e.g., 'charlie=T' we assume a priori
+    #  that they are illegal
+    #
+    dotargs <- list(...)
+    named <- if (is.null(names(dotargs))) 
+	           rep(FALSE, length(dotargs))
+             else (names(dotargs) != "")
+    if (any(named)) 
+        warning(paste("The following arguments to anova.coxphms(..)", 
+            "are invalid and dropped:", paste(deparse(dotargs[named]), 
+                collapse = ", ")))
+
+    allmod <- c(object, dotargs[!named])  # all the models
+    is.coxms <- sapply(dotargs, function(x) inherits(x, "coxphms"))
+    if (!all(is.coxms))
+        stop("All arguments must be multistate hazard models")    
+
+    ties <- sapply(allmod, function(x) x$method)
+    if (any(ties != ties[1]))
+        stop("all models must have the same ties option")
+stop("anova not yet available for multistate")
+    responses <- as.character(unlist(lapply(allmod, 
+				     function(x) deparse(formula(x)[[2]]))))
+    sameresp <- (responses == responses[1])
+    if (!all(sameresp)) {
+        allmod <- allmod[sameresp]
+        warning(paste("Models with response", deparse(responses[!sameresp]), 
+            "removed because response differs from", "model 1"))
+    }
+    nmodel <- length(allmod)
+    if (nmodel < 2) stop("must have more than one model")
+
+    # Check that they were all fit to the same data set.  This isn't perfect:
+    #  we could have distince data sets with the same number of rows in the
+    #  model frame (n), same number of unique id values, and same number of
+    #  events.  But it is unlikely.  The main reason for this is that one fit
+    #  had an extra variable that was missing on some subjects.
+    ns <- sapply(object, function(x) c(x$n, x$n.id, x$nevent))
+    if (!(is.matrix(ns))) # only occurs if a user messed with the object
+        stop("at least one model is missing the n, n.id, or nevent element")
+    if (any(apply(ns, 2, function(x) any(x != x[1])))) 
+        stop("models were not all fit to the same dataset")
+
+    # I can only handle models with the same stratification structure
+    stest <- sapply(allmod, function(x) identical(x$map, allmod[[1]]$smap))
+    if (!all(stest))
+        stop("not all models have the same structure of baseline hazards")
+    
+    # Models must be in increasing order of complexity
+    nvar <- sapply(allmod, function(x) length(coef(x)))
+    if (any(diff(nvar) < 1)) 
+        stop("models must be in increasing order of complexity")
+    # do the more complex nesting via variable names
+    for (i in 2:nmodel) {
+        indx <- match(rownames(allmod[i-1][["cmap"]]),
+                      rownames(allmod[i][["cmap"]]), nomatch=0)
+        if (any(indx==0))
+            stop(paste("model", i-1, "contains variables not in model", i))
+    }
+ 
+    # They all must have the same response
+    for (i in 2:nmodel)
+        if (!identical(allmod[[i]]$y, allmod[[i-1]]$y))
+            stop("all models must have the same response")
+
+    # Now for the real work
+    test <- df <- pval <- double(nmodel -1)
+    for (i in 2:nmodel)
+    
+    tfun <- function(x) paste(as.character(delete.response(terms(formula(x)))),
+                              collapse=' ')
+    variables <- lapply(object, tfun)
+    dimnames(table) <- list(1:nmodel, 
+			    c("loglik", "Chisq", "Df"))
+    title <- paste("Analysis of Deviance Table\n Cox model: response is ",
+		   responses[1]) 
+    topnote <- paste(" Model ", format(1:nmodel), ": ", variables, 
+		     sep = "", collapse = "\n")
+    if (!is.null(test)) {
+        table[['Pr(>|Chi|)']] <- pchisq(table$Chisq, table$Df, lower.tail=FALSE)
+        }
+    structure(table, heading = c(title, topnote), 
+			  class = c("anova", "data.frame"))
+} 
+

R/anova.coxphms.R

@@ -12,15 +12,25 @@ concordance.formula <- function(object, data,
     timewt <- match.arg(timewt)
     if (missing(ymin)) ymin <- NULL
     if (missing(ymax)) ymax <- NULL
-    
+    if (missing(object)) stop("a formula argument is required")
+    formula <- object  # clearer to read below
+
+    # make Surv(), strata() etc in a formula resolve to the survival namespace
+    newform <- removeDoubleColonSurv(formula)
+    if (!is.null(newform)) {
+        formula <- newform$formula
+        if (newform$newcall) Call$formula <- formula
+    }
+   
     index <- match(c("data", "weights", "subset", "na.action", 
                      "cluster"),
                    names(Call), nomatch=0)
     temp <- Call[c(1, index)]
     temp[[1L]] <-  quote(stats::model.frame)
     special <- c("strata", "cluster")
-    temp$formula <- if(missing(data)) terms(object, special)
-                    else              terms(object, special, data=data)
+    temp$formula <- if(missing(data)) terms(formula, special)
+                    else              terms(formula, special, data=data)
+
     mf <- eval(temp, parent.frame())  # model frame
     if (nrow(mf) ==0) stop("No (non-missing) observations")
     Terms <- terms(mf)

R/concordance.R

@@ -1,5 +1,4 @@
 # Automatically generated from the noweb directory
-#tt <- function(x) x
 coxph <- function(formula, data, weights, subset, na.action,
         init, control, ties= c("efron", "breslow", "exact"),
         singular.ok =TRUE,  robust,
@@ -9,6 +8,8 @@ coxph <- function(formula, data, weights, subset, na.action,
     missing.ties <- missing(ties) & missing(method) #see later multistate sect
     ties <- match.arg(ties)
     Call <- match.call()
+    if (missing(formula)) stop("a formula argument is required")
+    
     ## We want to pass any ... args to coxph.control, but not pass things
     ##  like "dats=mydata" where someone just made a typo.  The use of ...
     ##  is simply to allow things like "eps=1e6" with easier typing
@@ -20,14 +21,32 @@ coxph <- function(formula, data, weights, subset, na.action,
             stop(gettextf("Argument %s not matched", 
                           names(extraArgs)[indx==0L]), domain = NA)
     }
+    
+    # Gather any leftover arguments into a coxph.control call
+    # If there is a control argument, force a call to coxph.control to both
+    #  fill it out with all the elements and do sanity checks
     if (missing(control)) control <- coxph.control(...) 
+    else if (is.list(control)) control <- do.call(coxph.control, control)
+    else stop("control argument must be a list")
 
+    # make Surv(), strata() etc in a formula resolve to the survival namespace
+    if (is.list(formula)) {
+        newform <- removeDoubleColonSurv(formula[[1]])
+        if (!is.null(newform)) {
+            formula[[1]] <- newform$formula
+            if (newform$newcall) Call$formula <- newform$formula
+        }
+    } else {
+        newform <- removeDoubleColonSurv(formula)
+        if (!is.null(newform)) {
+            formula <- newform$formula
+            if (newform$newcall) Call$formula <- formula  #save the nicer version
+        }
+    }
     # Move any cluster() term out of the formula, and make it an argument
     #  instead.  This makes everything easier.  But, I can only do that with
     #  a local copy, doing otherwise messes up future use of update() on
     #  the model object for a user stuck in "+ cluster()" mode.
-    if (missing(formula)) stop("a formula argument is required")
-    
     ss <- "cluster"
     if (is.list(formula))
         Terms <- if (missing(data)) terms(formula[[1]], specials=ss) else
@@ -60,7 +79,6 @@ coxph <- function(formula, data, weights, subset, na.action,
     indx <- match(c("formula", "data", "weights", "subset", "na.action",
                     "cluster", "id", "istate"),
                   names(Call), nomatch=0) 
-    if (indx[1] ==0) stop("A formula argument is required")
     tform <- Call[c(1,indx)]  # only keep the arguments we wanted
     tform[[1L]] <- quote(stats::model.frame)  # change the function called
 
@@ -99,13 +117,6 @@ coxph <- function(formula, data, weights, subset, na.action,
     tform$formula <- if(missing(data)) terms(formula, special) else
                                       terms(formula, special, data=data)
 
-    # Make "tt" visible for coxph formulas, without making it visible elsewhere
-    if (!is.null(attr(tform$formula, "specials")$tt)) {
-        coxenv <- new.env(parent= environment(formula))
-        assign("tt", function(x) x, envir=coxenv)
-        environment(tform$formula) <- coxenv
-    }
-
     # okay, now evaluate the formula
     mf <- eval(tform, parent.frame())
     Terms <- terms(mf)
@@ -155,9 +166,14 @@ coxph <- function(formula, data, weights, subset, na.action,
             stop("multi-state models do not currently support pspline terms")
         if (length(attr(Terms, "specials")$ridge) >0)
             stop("multi-state models do not currently support ridge penalties")
-        if (!missing.ties) method <- ties <- "breslow"
+        if (missing.ties) method <- ties <- "breslow"
     }
 
+    # the code was never designed for multiple fraily terms, but of course
+    #  someone tried it
+    if (length(attr(Terms, "specials")$frailty) >1)
+            stop("multiple frailty terms are not supported")
+
     if (control$timefix) Y <- aeqSurv(Y)
     if (length(attr(Terms, 'variables')) > 2) { # a ~1 formula has length 2
         ytemp <- innerterms(formula[1:2])
@@ -191,10 +207,6 @@ coxph <- function(formula, data, weights, subset, na.action,
         if (!hasinteractions) dropterms <- stemp$terms 
     } else istrat <- NULL
 
-    if (hasinteractions && multi)
-        stop("multi-state coxph does not support strata*covariate interactions")
-
-
     timetrans <- attr(Terms, "specials")$tt
     if (missing(tt)) tt <- NULL
     if (length(timetrans)) {

R/coxph.R

@@ -6,12 +6,18 @@ coxph.control <- function(eps=1e-9,
 			  iter.max=20,
 			  toler.inf= sqrt(eps), outer.max=10,
                           timefix =TRUE) {
-    if (iter.max <0) stop("Invalid value for iterations")
-    if (eps <=0) stop ("Invalid convergence criteria")
+    if (!is.numeric(iter.max) ||iter.max <0) stop("Invalid value for iterations")
+    if (!is.numeric(eps) || eps <=0) stop ("Invalid convergence criteria")
+    if (!is.numeric(toler.chol) || toler.chol <=0)
+        stop("invalid value for toler.chol")
+    if (!is.numeric(eps) || eps <=0) stop("eps must be > 0")
     if (eps <= toler.chol) 
 	    warning("For numerical accuracy, tolerance should be < eps")
-    if (toler.inf <=0) stop ("The inf.warn setting must be >0")
+    if (!is.numeric(toler.inf) || toler.inf <=0) 
+        stop ("The toler.inf setting must be >0")
     if (!is.logical(timefix)) stop("timefix must be TRUE or FALSE")
+    if (!is.numeric(outer.max) || outer.max <=0) 
+        stop("invalid value for outer.max")
     list(eps=eps, toler.chol=toler.chol, iter.max=as.integer(iter.max), 
 	 toler.inf=toler.inf, outer.max=as.integer(outer.max), 
          timefix=timefix)