1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
|
# Automatically generated from all.nw using noweb
survfit <- function(formula, ...) {
UseMethod("survfit", formula)
}
dim.survfit <- function(x) {
if (is.null(x$strata)) {
if (is.matrix(x$surv)) ncol(x$surv)
else 1
}
else {
nr <- length(x$strata)
if (is.matrix(x$xurv)) c(nr, ncol(x$surv))
else nr
}
}
"[.survfit" <- function(x, ..., drop=TRUE) {
nmatch <- function(indx, target) {
# This function lets R worry about character, negative, or logical subscripts
# It always returns a set of positive integer indices
temp <- 1:length(target)
names(temp) <- target
temp[indx]
}
if (missing(..1)) i<- NULL else i <- ..1
if (missing(..2)) j<- NULL else j <- ..2
if (is.null(i) && is.null(j)) return (x) #no subscripts present!
if (!is.matrix(x$surv) && !is.null(j))
stop("survfit object does not have 2 dimensions")
if (is.null(x$strata)) {
if (is.matrix(x$surv)) {
if (is.null(j) && !is.null(i)) j <- i #special case noted above
x$surv <- x$surv[,j,drop=drop]
if (!is.null(x$std.err)) x$std.err <- x$std.err[,j,drop=drop]
if (!is.null(x$upper)) x$upper <- x$upper[,j,drop=drop]
if (!is.null(x$lower)) x$lower <- x$lower[,j,drop=drop]
if (!is.null(x$cumhaz)) x$cumhaz <- x$cumhaz[,j,drop=drop]
}
else warning("survfit object has only a single survival curve")
}
else {
if (is.null(i)) keep <- seq(along.with=x$time)
else {
indx <- nmatch(i, names(x$strata)) #strata to keep
if (any(is.na(indx)))
stop(paste("strata",
paste(i[is.na(indx)], collapse=' '),
'not matched'))
# Now, indx may not be in order: some can use curve[3:2] to reorder
# The list/unlist construct will reorder the data
temp <- rep(1:length(x$strata), x$strata)
keep <- unlist(lapply(indx, function(x) which(temp==x)))
if (length(indx) <=1 && drop) x$strata <- NULL
else x$strata <- x$strata[i]
x$n <- x$n[indx]
x$time <- x$time[keep]
x$n.risk <- x$n.risk[keep]
x$n.event <- x$n.event[keep]
x$n.censor<- x$n.censor[keep]
if (!is.null(x$enter)) x$enter <- x$enter[keep]
}
if (is.matrix(x$surv)) {
# If the curve has been selected by strata and keep has only
# one row, we don't want to lose the second subscript too
if (!is.null(i) && (is.null(j) ||length(j) >1)) drop <- FALSE
if (is.null(j)) {
x$surv <- x$surv[keep,,drop=drop]
if (!is.null(x$std.err))
x$std.err <- x$std.err[keep,,drop=drop]
if (!is.null(x$upper)) x$upper <-x$upper[keep,,drop=drop]
if (!is.null(x$lower)) x$lower <-x$lower[keep,,drop=drop]
if (!is.null(x$cumhaz)) x$cumhaz <-x$cumhaz[keep,,drop=drop]
}
else {
x$surv <- x$surv[keep,j, drop=drop]
if (!is.null(x$std.err))
x$std.err <- x$std.err[keep,j, drop=drop]
if (!is.null(x$upper)) x$upper <- x$upper[keep,j, drop=drop]
if (!is.null(x$lower)) x$lower <- x$lower[keep,j, drop=drop]
if (!is.null(x$cumhaz)) x$cumhaz <- x$cumhaz[keep,j, drop=drop]
}
}
else {
x$surv <- x$surv[keep]
if (!is.null(x$std.err)) x$std.err <- x$std.err[keep]
if (!is.null(x$upper)) x$upper <- x$upper[keep]
if (!is.null(x$lower)) x$lower <- x$lower[keep]
if (!is.null(x$cumhaz)) x$cumhaz <- x$cumhaz[keep]
}
}
x
}
survfit.formula <- function(formula, data, weights, subset,
na.action, etype, id, istate, ...) {
Call <- match.call()
Call[[1]] <- as.name('survfit') #make nicer printout for the user
# create a copy of the call that has only the arguments we want,
# and use it to call model.frame()
indx <- match(c('formula', 'data', 'weights', 'subset','na.action',
'istate', 'id', "etype"), names(Call), nomatch=0)
#It's very hard to get the next error message other than malice
# eg survfit(wt=Surv(time, status) ~1)
if (indx[1]==0) stop("a formula argument is required")
temp <- Call[c(1, indx)]
temp[[1]] <- as.name("model.frame")
m <- eval.parent(temp)
Terms <- terms(formula, c("strata", "cluster"))
ord <- attr(Terms, 'order')
if (length(ord) & any(ord !=1))
stop("Interaction terms are not valid for this function")
n <- nrow(m)
Y <- model.extract(m, 'response')
if (!is.Surv(Y)) stop("Response must be a survival object")
casewt <- model.extract(m, "weights")
if (is.null(casewt)) casewt <- rep(1,n)
if (!is.null(attr(Terms, 'offset'))) warning("Offset term ignored")
id <- model.extract(m, 'id')
istate <- model.extract(m,"istate")
temp <- untangle.specials(Terms, "cluster")
if (length(temp$vars)>0) {
if (length(temp$vars) > 1) stop("can not have two cluster terms")
if (!is.null(id)) stop("can not have both a cluster term and an id variable")
id <- m[[temp$vars]]
Terms <- Terms[-temp$terms]
}
ll <- attr(Terms, 'term.labels')
if (length(ll) == 0) X <- factor(rep(1,n)) # ~1 on the right
else X <- strata(m[ll])
if (!is.Surv(Y)) stop("y must be a Surv object")
# Backwards support for the now-depreciated etype argument
etype <- model.extract(m, "etype")
if (!is.null(etype)) {
if (attr(Y, "type") == "mcounting" ||
attr(Y, "type") == "mright")
stop("cannot use both the etype argument and mstate survival type")
if (length(istate))
stop("cannot use both the etype and istate arguments")
status <- Y[,ncol(Y)]
etype <- as.factor(etype)
temp <- table(etype, status==0)
if (all(rowSums(temp==0) ==1)) {
# The user had a unique level of etype for the censors
newlev <- levels(etype)[order(-temp[,2])] #censors first
}
else newlev <- c(" ", levels(etype)[temp[,1] >0])
status <- factor(ifelse(status==0,0, as.numeric(etype)),
labels=newlev)
if (attr(Y, 'type') == "right")
Y <- Surv(Y[,1], status, type="mstate")
else if (attr(Y, "type") == "counting")
Y <- Surv(Y[,1], Y[,2], status, type="mstate")
else stop("etype argument incompatable with survival type")
}
# At one point there were lines here to round the survival
# times to a certain number of digits. This approach worked
# almost all the time, but only almost. The better logic is
# now in the individual compuation routines
if (attr(Y, 'type') == 'left' || attr(Y, 'type') == 'interval')
temp <- survfitTurnbull(X, Y, casewt, ...)
else if (attr(Y, 'type') == "right" || attr(Y, 'type')== "counting")
temp <- survfitKM(X, Y, casewt, ...)
else if (attr(Y, 'type') == "mright" || attr(Y, "type")== "mcounting")
temp <- survfitCI(X, Y, weights=casewt, id=id, istate=istate, ...)
else {
# This should never happen
stop("unrecognized survival type")
}
if (is.null(temp$states)) class(temp) <- 'survfit'
else class(temp) <- c("survfitms", "survfit")
if (!is.null(attr(m, 'na.action')))
temp$na.action <- attr(m, 'na.action')
temp$call <- Call
temp
}
|