diff --git a/R/survfit.coxphms.R b/R/survfit.coxphms.R index bc0b5adf..aece724f 100644 --- a/R/survfit.coxphms.R +++ b/R/survfit.coxphms.R @@ -1,7 +1,7 @@ # Automatically generated from the noweb directory survfit.coxphms <- function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, - stype=2, ctype, + stype=2, ctype, conf.type=c("log", "log-log", "plain", "none", "logit", "arcsin"), censor=TRUE, start.time, id, influence=FALSE, na.action=na.pass, type, p0=NULL, time0=FALSE, ...) { @@ -12,9 +12,9 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, se.fit <- FALSE #still to do if (missing(newdata)) stop("multi-state survival requires a newdata argument") - if (!missing(id)) + if (!missing(id)) stop("using a covariate path is not supported for multi-state") - temp <- object$smap["(Baseline)",] + temp <- object$smap["(Baseline)",] baselinecoef <- rbind(temp, coef= 1.0) phbase <- rep(FALSE, nrow(object$cmap)) if (any(duplicated(temp))) { @@ -24,7 +24,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, # There might not be such rows, by the way. pattern <- "^ph\\([0-9]+:[0-9]+\\)$" cname <- rownames(object$cmap) - phbase <- grepl(pattern, cname) # this row points to a "ph" coef + phbase <- grepl(pattern, cname) # this row points to a "ph" coef for (i in which(phbase)) { # Say that this row (i) of cmap had label ph(1:4), and contains # elements 0,0,0,0,0, 8,9. @@ -36,7 +36,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, baselinecoef[2, j>0] <- exp(object$coef[j]) } } - + # process options, set up Y and the model frame for the original data Terms <- terms(object) robust <- !is.null(object$naive.var) # did the coxph model use robust var? @@ -163,7 +163,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, istate <- istate[-toss] } } - + # expansion of the X matrix with stacker, set up shared hazards # Rebuild istate using the survcheck routine, as a double check # that the data set hasn't been modified @@ -175,7 +175,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, else { # if istate has unused levels, mcheck$istate won't have them so they # need to be dropped. - istate <- factor(istate, object$states) + istate <- factor(istate, object$states) # a new level in state should only happen if someone has mucked up the # data set used in the coxph fit if (any(is.na(istate))) stop("unrecognized initial state, data changed?") @@ -190,14 +190,14 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, if (is.null(strata)) tempstrat <- rep(1L, nrow(Y)) else tempstrat <- strata - cifit <- survfitAJ(as.factor(tempstrat), Y, weights, - id= oldid, istate = istate, se.fit=FALSE, + cifit <- survfitAJ(as.factor(tempstrat), Y, weights, + id= oldid, istate = istate, se.fit=FALSE, start.time=start.time, p0=p0, time0= time0) # For computing the actual estimates it is easier to work with an # expanded data set. # Replicate actions found in the coxph-multi-X chunk - # Note the dropzero=FALSE argument: if there is a transition with no + # Note the dropzero=FALSE argument: if there is a transition with no # covariates we still need it expanded; this differs from coxph. # A second differnence is tstrata: force stacker to think that every # transition is a unique hazard, so that it does proper expansion. @@ -378,7 +378,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, # x2 will have one row per desired curve and one col per 'normal' covariate. risk2 <- exp(x2 %*% ifelse(is.na(temp), 0, temp) - xcenter) # risk2 has a risk score with rows= curve and cols= transition - # make the expansion map. + # make the expansion map. # The H matrices we will need are nstate by nstate, at each time, with # elements that are non-zero only for observed transtions. states <- object$states @@ -393,7 +393,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, ny <- ncol(Y) if (is.null(strata)) { fit <- multihaz(Y, X, position, weights, risk, istrat, ctype, stype, - baselinecoef, hfill, x2, risk2, varmat, nstate, se.fit, + baselinecoef, hfill, x2, risk2, varmat, nstate, se.fit, cifit$p0, cifit$time) cifit$pstate <- fit$pstate cifit$cumhaz <- fit$cumhaz @@ -432,8 +432,8 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, class(cifit) <- c("survfitcoxms", "survfitms", "survfit") cifit } -# Compute the hazard and survival functions -multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, +# Compute the hazard and survival functions +multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, bcoef, hfill, x2, risk2, vmat, nstate, se.fit, p0, utime) { ny <- ncol(y) sort2 <- order(istrat, y[,ny-1L]) -1L @@ -443,22 +443,22 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, # this returns all of the counts we might desire. if (ny ==2) { fit <- .Call(Ccoxsurv1, utime, y, weight, sort2, istrat, x, risk) - cn <- fit$count - dim(cn) <- c(length(utime), fit$ntrans, 10) + cn <- fit$count + dim(cn) <- c(length(utime), fit$ntrans, 10) } - else { + else { sort1 <- order(istrat, y[,1]) -1L - fit <- .Call(Ccoxsurv2, utime, y, weight, sort1, sort2, position, + fit <- .Call(Ccoxsurv2, utime, y, weight, sort1, sort2, position, istrat, x, risk) - cn <- fit$count - dim(cn) <- c(length(utime), fit$ntrans, 12) + cn <- fit$count + dim(cn) <- c(length(utime), fit$ntrans, 12) } # cn is returned as a matrix since there is an allocMatrix C macro, but # no allocArray macro. So we first reset the dimensions. # The first dimension is time # Second is the transition, same order as columns of bcoef # Third is the count type: 1-3 = at risk (unweighted, with case weights, - # with casewt * risk wt), 4-6 = events (unweighted, case, risk), + # with casewt * risk wt), 4-6 = events (unweighted, case, risk), # 7-8 = censored events, 9-10 = censored, 11-12 = Efron # We will use events/(at risk) = cn[,,5]/cn[,,3] a few lines below; avoid 0/0 @@ -471,8 +471,8 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, denom1 <- ifelse(none.atrisk, 1, cn[,,3]) # avoid a later 0/0 denom2 <- ifelse(none.atrisk, 1, cn[,,3]^2) } else { - denom1 <- ifelse(none.atrisk, 1, cn[,,9]) - denom2 <- ifelse(none.atrisk, 1, cn[,,10]) + denom1 <- ifelse(none.atrisk, 1, 1/cn[,,9]) + denom2 <- ifelse(none.atrisk, 1, 1/cn[,,10]) } # We want to avoid 0/0. If there is no one at risk (denominator) then @@ -488,10 +488,10 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, else atrisk <- cn[,,9] %*% design basehaz <- events/ifelse(atrisk<=0, 1, atrisk) hazard <- basehaz[,bcoef[1,]] * rep(bcoef[2,], each=nrow(basehaz)) - } + } else { if (ctype==1) hazard <- cn[,,5]/ifelse(cn[,,3]<=0, 1, cn[,,3]) - else hazard <- cn[,,5]/ifelse(cn[,,9] <=0, 1, cn[,,9]) + else hazard <- cn[,,5] * ifelse(cn[,,9] <=0, 1, cn[,,9]) } # Expand the result, one "hazard set" for each row of x2 @@ -499,7 +499,7 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, h2 <- array(0, dim=c(nrow(hazard), nx2, ncol(hazard))) S <- double(nstate) # survival at the current time S2 <- array(0, dim=c(nrow(hazard), nx2, nstate)) - + H <- matrix(0, nstate, nstate) if (stype==2) { H[hfill] <- colMeans(hazard) # dummy H to drive esetup diff --git a/noweb/code.pdf b/noweb/code.pdf new file mode 100644 index 00000000..3698b1e6 Binary files /dev/null and b/noweb/code.pdf differ diff --git a/noweb/coxsurv3.Rnw b/noweb/coxsurv3.Rnw index 28736ff1..4e5fb932 100644 --- a/noweb/coxsurv3.Rnw +++ b/noweb/coxsurv3.Rnw @@ -5,7 +5,7 @@ particularly the variance. <>= survfit.coxphms <- function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, - stype=2, ctype, + stype=2, ctype, conf.type=c("log", "log-log", "plain", "none", "logit", "arcsin"), censor=TRUE, start.time, id, influence=FALSE, na.action=na.pass, type, p0=NULL, time0=FALSE, ...) { @@ -16,9 +16,9 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, se.fit <- FALSE #still to do if (missing(newdata)) stop("multi-state survival requires a newdata argument") - if (!missing(id)) + if (!missing(id)) stop("using a covariate path is not supported for multi-state") - temp <- object$smap["(Baseline)",] + temp <- object$smap["(Baseline)",] baselinecoef <- rbind(temp, coef= 1.0) phbase <- rep(FALSE, nrow(object$cmap)) if (any(duplicated(temp))) { @@ -28,7 +28,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, # There might not be such rows, by the way. pattern <- "^ph\\([0-9]+:[0-9]+\\)$" cname <- rownames(object$cmap) - phbase <- grepl(pattern, cname) # this row points to a "ph" coef + phbase <- grepl(pattern, cname) # this row points to a "ph" coef for (i in which(phbase)) { # Say that this row (i) of cmap had label ph(1:4), and contains # elements 0,0,0,0,0, 8,9. @@ -40,7 +40,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, baselinecoef[2, j>0] <- exp(object$coef[j]) } } - + # process options, set up Y and the model frame for the original data <> <> @@ -62,7 +62,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, istate <- istate[-toss] } } - + # expansion of the X matrix with stacker, set up shared hazards <> @@ -77,7 +77,7 @@ function(formula, newdata, se.fit=FALSE, conf.int=.95, individual=FALSE, class(cifit) <- c("survfitcoxms", "survfitms", "survfit") cifit } -@ +@ The third line \code{as.name('survfit')} causes the printout to say `survfit' instead of `survfit.coxph'. %' @@ -98,7 +98,7 @@ if (is.null(istate)) istate <- mcheck$istate else { # if istate has unused levels, mcheck$istate won't have them so they # need to be dropped. - istate <- factor(istate, object$states) + istate <- factor(istate, object$states) # a new level in state should only happen if someone has mucked up the # data set used in the coxph fit if (any(is.na(istate))) stop("unrecognized initial state, data changed?") @@ -113,14 +113,14 @@ if (is.null(weights)) weights <- rep(1.0, nrow(Y)) if (is.null(strata)) tempstrat <- rep(1L, nrow(Y)) else tempstrat <- strata -cifit <- survfitAJ(as.factor(tempstrat), Y, weights, - id= oldid, istate = istate, se.fit=FALSE, +cifit <- survfitAJ(as.factor(tempstrat), Y, weights, + id= oldid, istate = istate, se.fit=FALSE, start.time=start.time, p0=p0, time0= time0) # For computing the actual estimates it is easier to work with an # expanded data set. # Replicate actions found in the coxph-multi-X chunk -# Note the dropzero=FALSE argument: if there is a transition with no +# Note the dropzero=FALSE argument: if there is a transition with no # covariates we still need it expanded; this differs from coxph. # A second differnence is tstrata: force stacker to think that every # transition is a unique hazard, so that it does proper expansion. @@ -146,7 +146,7 @@ if (robust & length(cluster)==0) cluster <- oldid Fix up the X matrix to avoid huge values. In the single state case this is fairly straightforward: use $(X-1m')\beta = X\beta - m'\beta$ where -$m$ is the vector of centering constants found in the +$m$ is the vector of centering constants found in the \code{object\$means} component. However, in multi-state there will often be covariates that are part of one transition but not another, and if one of them is wild we will want different @@ -173,14 +173,14 @@ else { keep <- !grepl("frailty(", dimnames(X)[[2]], fixed=TRUE) X <- X[,keep, drop=F] } - + if (is.null(offset)) risk <- c(exp(X%*% beta - xcenter)) else risk <- c(exp(X%*% beta + offset - xcenter)) } -@ +@ The survfit.coxph-setup3 chunk, shared with single state Cox models, has created -an mf2 model frame and an x2 matrix. +an mf2 model frame and an x2 matrix. For multi-state, we ignore any strata variables in mf2. Create a matrix of risk scores, number of subjects by number of transitions. Different transitions often have different coefficients, so there is a risk @@ -198,7 +198,7 @@ temp <- coef(object, matrix=TRUE)[!phbase,,drop=FALSE] # ignore missing coefs # x2 will have one row per desired curve and one col per 'normal' covariate. risk2 <- exp(x2 %*% ifelse(is.na(temp), 0, temp) - xcenter) # risk2 has a risk score with rows= curve and cols= transition -@ +@ At this point we have several parts to keep straight. The data set has been expanded into a new X and Y. @@ -222,7 +222,7 @@ The coxph routine fits all strata and transitions at once, since the loglik is a sum over strata. This routine does each stratum separately. <>= -# make the expansion map. +# make the expansion map. # The H matrices we will need are nstate by nstate, at each time, with # elements that are non-zero only for observed transtions. states <- object$states @@ -237,7 +237,7 @@ if (individual) { ny <- ncol(Y) if (is.null(strata)) { fit <- multihaz(Y, X, position, weights, risk, istrat, ctype, stype, - baselinecoef, hfill, x2, risk2, varmat, nstate, se.fit, + baselinecoef, hfill, x2, risk2, varmat, nstate, se.fit, cifit$p0, cifit$time) cifit$pstate <- fit$pstate cifit$cumhaz <- fit$cumhaz @@ -275,19 +275,19 @@ cifit$newdata <- newdata Finally, a routine that does all the actual work. \begin{itemize} - \item The first 5 variables are for the data set that the Cox model was built - on: y, x, position, risk score, istrat. + \item The first 5 variables are for the data set that the Cox model was built + on: y, x, position, risk score, istrat. Position is a flag for each obs. Is it the first of a connected string - such as (10, 12) (12,19) (19,21), the last of such a string, both, + such as (10, 12) (12,19) (19,21), the last of such a string, both, or neither. 1*first + 2*last. This affects whether an obs is labeled as censored or not in user printout, nothing else. (That part has actually already been done via the survfitAJ call.) - \item x2 and risk2 are the covariates and risk scores for the predicted + \item x2 and risk2 are the covariates and risk scores for the predicted values. These do not involve any ph(a:b) coefficients. \item baselinecoef encodes shared hazards - \item hfill control mapping from fitted hazards to + \item hfill control mapping from fitted hazards to transitions and probabilities - \item p0 will be NULL if the user did not specifiy it. + \item p0 will be NULL if the user did not specifiy it. \item vmat is only needed for standard errors \item utime is the set of time points desired \end{itemize} @@ -301,8 +301,8 @@ and the 3 individual hazards are 1*base, 1.3*base and .4*base. If there are no shared hazards this can be computed more simply of course. <>= -# Compute the hazard and survival functions -multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, +# Compute the hazard and survival functions +multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, bcoef, hfill, x2, risk2, vmat, nstate, se.fit, p0, utime) { ny <- ncol(y) sort2 <- order(istrat, y[,ny-1L]) -1L @@ -312,22 +312,22 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, # this returns all of the counts we might desire. if (ny ==2) { fit <- .Call(Ccoxsurv1, utime, y, weight, sort2, istrat, x, risk) - cn <- fit$count - dim(cn) <- c(length(utime), fit$ntrans, 10) + cn <- fit$count + dim(cn) <- c(length(utime), fit$ntrans, 10) } - else { + else { sort1 <- order(istrat, y[,1]) -1L - fit <- .Call(Ccoxsurv2, utime, y, weight, sort1, sort2, position, + fit <- .Call(Ccoxsurv2, utime, y, weight, sort1, sort2, position, istrat, x, risk) - cn <- fit$count - dim(cn) <- c(length(utime), fit$ntrans, 12) + cn <- fit$count + dim(cn) <- c(length(utime), fit$ntrans, 12) } # cn is returned as a matrix since there is an allocMatrix C macro, but # no allocArray macro. So we first reset the dimensions. # The first dimension is time # Second is the transition, same order as columns of bcoef # Third is the count type: 1-3 = at risk (unweighted, with case weights, - # with casewt * risk wt), 4-6 = events (unweighted, case, risk), + # with casewt * risk wt), 4-6 = events (unweighted, case, risk), # 7-8 = censored events, 9-10 = censored, 11-12 = Efron # We will use events/(at risk) = cn[,,5]/cn[,,3] a few lines below; avoid 0/0 @@ -340,8 +340,8 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, denom1 <- ifelse(none.atrisk, 1, cn[,,3]) # avoid a later 0/0 denom2 <- ifelse(none.atrisk, 1, cn[,,3]^2) } else { - denom1 <- ifelse(none.atrisk, 1, cn[,,9]) - denom2 <- ifelse(none.atrisk, 1, cn[,,10]) + denom1 <- ifelse(none.atrisk, 1, 1/cn[,,9]) + denom2 <- ifelse(none.atrisk, 1, 1/cn[,,10]) } # We want to avoid 0/0. If there is no one at risk (denominator) then @@ -357,10 +357,10 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, else atrisk <- cn[,,9] %*% design basehaz <- events/ifelse(atrisk<=0, 1, atrisk) hazard <- basehaz[,bcoef[1,]] * rep(bcoef[2,], each=nrow(basehaz)) - } + } else { if (ctype==1) hazard <- cn[,,5]/ifelse(cn[,,3]<=0, 1, cn[,,3]) - else hazard <- cn[,,5]/ifelse(cn[,,9] <=0, 1, cn[,,9]) + else hazard <- cn[,,5] * ifelse(cn[,,9] <=0, 1, cn[,,9]) } # Expand the result, one "hazard set" for each row of x2 @@ -368,7 +368,7 @@ multihaz <- function(y, x, position, weight, risk, istrat, ctype, stype, h2 <- array(0, dim=c(nrow(hazard), nx2, ncol(hazard))) S <- double(nstate) # survival at the current time S2 <- array(0, dim=c(nrow(hazard), nx2, nstate)) - + H <- matrix(0, nstate, nstate) if (stype==2) { H[hfill] <- colMeans(hazard) # dummy H to drive esetup diff --git a/src/coxsurv1.c b/src/coxsurv1.c index 13454fd4..19c5dbfc 100644 --- a/src/coxsurv1.c +++ b/src/coxsurv1.c @@ -18,15 +18,15 @@ ** routine will be called separatedly for individual strata, with otime the ** event times in that stratum. Ditto for a multistate model that had ** external strata. -** +** ** otime: vector of output times. All the transitions will get reports at -** these time points. This fcn is called for all of the +** these time points. This fcn is called for all of the ** transitions at once, sorted by transition, ** but called separately for any strata() groups. ** y : survival response, two column ** weight: observation weight ** sort2: sort index for the survival time -** trans: the data set is stacked: all the data for transition 1, then +** trans: the data set is stacked: all the data for transition 1, then ** transition 2, etc for a multi-state model (th ** xmat2: covariates ** risk2: risk score @@ -35,7 +35,7 @@ ** matrices with xbar for those at risk, and the sum of x for ** terminal events at the current time. ** -** For the weighted counts, number at risk != entries - exits. Someone with +** For the weighted counts, number at risk != entries - exits. Someone with ** a sequence of (1,2)(2,5)(5,6) will have 1 entry and 1 exit, but they might ** have 3 changes of risk score due to time-dependent covariates. ** n0-3 has to count all the changes, while n8-n9 (only used in printout) @@ -54,9 +54,9 @@ #include "survproto.h" #include -SEXP coxsurv1(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort22, +SEXP coxsurv1(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort22, SEXP trans2, SEXP xmat2, SEXP risk2) { - + int i, i2, k, person2, itrans; int nused, ntrans, ntime, irow, ii, jj; double *tstop, *status, *wt, *otime; @@ -67,8 +67,8 @@ SEXP coxsurv1(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort22, int nvar; /* number of covariates */ double *xsum1, /* a weighted sum, for computing xbar */ *xsum2; - - static const char *outnames[]={"ntrans", "count", + + static const char *outnames[]={"ntrans", "count", "xbar", "xsum2", ""}; SEXP rlist; double n[12]; @@ -100,29 +100,29 @@ SEXP coxsurv1(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort22, if (trans[i2] != itrans) { ntrans++; itrans = trans[i2]; - } - } - + } + } + /* Allocate memory for the working matrices. */ xsum1 = (double *) ALLOC(2*nvar, sizeof(double)); xsum2 = xsum1 + nvar; - + /* Allocate memory for returned objects: ntime*ntrans copies of n, xsum1, and xsum2 */ PROTECT(rlist = mkNamed(VECSXP, outnames)); irow = ntime*ntrans; rstrat = REAL(SET_VECTOR_ELT(rlist, 0, allocVector(REALSXP, 1))); - rn = dmatrix(REAL(SET_VECTOR_ELT(rlist, 1, + rn = dmatrix(REAL(SET_VECTOR_ELT(rlist, 1, allocMatrix(REALSXP, irow, 10))), irow, 10); - rx1 = dmatrix(REAL(SET_VECTOR_ELT(rlist, 2, + rx1 = dmatrix(REAL(SET_VECTOR_ELT(rlist, 2, allocMatrix(REALSXP, irow, nvar))), irow, nvar); - rx2 = dmatrix(REAL(SET_VECTOR_ELT(rlist, 3, + rx2 = dmatrix(REAL(SET_VECTOR_ELT(rlist, 3, allocMatrix(REALSXP, irow, nvar))), irow, nvar); - + R_CheckUserInterrupt(); /*check for control-C */ - /* now add up all the sums + /* now add up all the sums ** All this is done backwards in time. The logic is a bit easier, and ** the computation is numerically more stable (fewer subtractions). ** One by one for the desired output times "otime". @@ -130,15 +130,15 @@ SEXP coxsurv1(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort22, ** set when we cross their ending time. Also add them to the "censored" ** count. Don't add any who will be removed before otime to either ** count, however. - ** - ** 2. While tstop==otime, add them to the risk set, and count the + ** + ** 2. While tstop==otime, add them to the risk set, and count the ** observation with repect to n3- n7. ** ** In the code position2 is the current index into the sort2 vector, ** and i2 is the current value of sort2. */ rstrat[0] = ntrans; /* single element, number of transitions found */ - person2= nused-1; + person2= nused-1; irow = (ntime*ntrans); /* row of output objects */ for (ii =0; ii=0; jj--) { /* one by one through the times */ dtime = otime[jj]; - for (k=3; k<8; k++) n[k]=0; /* counts are only for this interval*/ + for (k=3; k<10; k++) n[k]=0; /* counts are only for this interval*/ /* Step 1 */ for(; person2 >=0 && trans[person2]==itrans; person2--) { @@ -156,7 +156,7 @@ SEXP coxsurv1(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort22, n[0]++; n[1] += wt[i2]; n[2] += wt[i2] * risk[i2]; - for (k=0; k -SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, +SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, SEXP sindex2, SEXP trans2, SEXP xmat2, SEXP risk2) { - + int i, i1, i2, k, person1, person2, itrans; int nused, ntrans, ntime, irow, ii, jj; double *tstart=0, *tstop, *status, *wt, *otime; @@ -60,10 +60,10 @@ SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, int nvar; /* number of covariates */ double *xsum1, /* a weighted sum, for computing xbar */ *xsum2; - + int *atrisk; - static const char *outnames[]={"ntrans", "count", + static const char *outnames[]={"ntrans", "count", "xbar", "xsum2", ""}; SEXP rlist; double n[12]; @@ -99,31 +99,31 @@ SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, if (trans[i2] != itrans) { ntrans++; itrans = trans[i2]; - } - } - + } + } + /* Allocate memory for the working matrices. */ xsum1 = (double *) ALLOC(2*nvar, sizeof(double)); xsum2 = xsum1 + nvar; atrisk = (int *) ALLOC(nused, sizeof(int)); for (i=0; i otime, remove any obs currently at risk from n0-n2. @@ -142,7 +142,7 @@ SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, ** and sort2 vectors, and i1/i2 are the current values of those vectors. */ rstrat[0] = ntrans; /* single element, number of transitions found */ - person1 = nused-1; person2= nused-1; /* person1 tracks start times */ + person1 = nused-1; person2= nused-1; /* person1 tracks start times */ irow = (ntime*ntrans); /* row of output objects */ for (ii =0; ii1 && status[i2]==0) { @@ -194,7 +194,7 @@ SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, i1 = sort1[person1]; if (tstart[i1] < dtime) break; if (atrisk[i1]) { /* remove them from risk set */ - n[0]--; + n[0]--; if (n[0] ==0) { n[1] =0; n[2] =0; @@ -203,7 +203,7 @@ SEXP coxsurv2(SEXP otime2, SEXP y2, SEXP weight2, SEXP sort12, SEXP sort22, else { n[1] -= wt[i1]; n[2] -= wt[i1]*risk[i1]; - for (k=0; k