diff -Nru survival-3.1-8/debian/changelog survival-3.1-11/debian/changelog --- survival-3.1-8/debian/changelog 2020-03-13 13:52:36.000000000 +0000 +++ survival-3.1-11/debian/changelog 2020-03-13 13:52:37.000000000 +0000 @@ -1,14 +1,21 @@ -survival (3.1-8-1cran1xenial0) xenial; urgency=medium +survival (3.1-11-1cran1xenial0) xenial; urgency=medium * Compilation for Ubuntu 16.04.6 LTS - -- Michael Rutter Sat, 14 Dec 2019 15:33:31 +0000 + -- Michael Rutter Fri, 13 Mar 2020 13:50:11 +0000 + +survival (3.1-11-1cran1) testing; urgency=low + + * cran2deb svn: 362M with DB version 1. + + -- cran2deb4ubuntu Sat, 07 Mar 2020 13:35:02 -0500 + survival (3.1-8-1cran1) testing; urgency=low * cran2deb svn: 362M with DB version 1. - -- cran2deb4ubuntu Sat, 07 Dec 2019 08:56:36 -0500 + -- cran2deb4ubuntu Sat, 07 Dec 2019 08:56:48 -0500 survival (3.1-7-1cran1) testing; urgency=low diff -Nru survival-3.1-8/debian/control survival-3.1-11/debian/control --- survival-3.1-8/debian/control 2020-03-13 13:52:36.000000000 +0000 +++ survival-3.1-11/debian/control 2020-03-13 13:52:37.000000000 +0000 @@ -17,6 +17,7 @@ Cox models, and parametric accelerated failure time models. . Author: Terry M Therneau [aut, cre], Thomas Lumley [ctb, trl] (original - S->R port and R maintainer until 2009) + S->R port and R maintainer until 2009), Atkinson Elizabeth [ctb], + Crowson Cynthia [ctb] . Maintainer: Terry M Therneau diff -Nru survival-3.1-8/debian/copyright survival-3.1-11/debian/copyright --- survival-3.1-8/debian/copyright 2020-03-13 13:52:36.000000000 +0000 +++ survival-3.1-11/debian/copyright 2020-03-13 13:52:37.000000000 +0000 @@ -2,9 +2,10 @@ automatically using cran2deb4ubuntu by cran2deb4ubuntu . -The original GNU R package is Copyright (C) 2019 Terry M Therneau [aut, +The original GNU R package is Copyright (C) 2020 Terry M Therneau [aut, cre], Thomas Lumley [ctb, trl] (original S->R port and R maintainer -until 2009) and possibly others. +until 2009), Atkinson Elizabeth [ctb], Crowson Cynthia [ctb] and +possibly others. The original GNU R package is maintained by Terry M Therneau and was obtained from: diff -Nru survival-3.1-8/DESCRIPTION survival-3.1-11/DESCRIPTION --- survival-3.1-8/DESCRIPTION 2019-12-03 17:30:02.000000000 +0000 +++ survival-3.1-11/DESCRIPTION 2020-03-07 12:10:03.000000000 +0000 @@ -2,8 +2,7 @@ Maintainer: Terry M Therneau Priority: recommended Package: survival -Version: 3.1-8 -Date: 2019-11-18 +Version: 3.1-11 Depends: R (>= 3.4.0) Imports: graphics, Matrix, methods, splines, stats, utils LazyData: Yes @@ -13,7 +12,9 @@ email="therneau.terry@mayo.edu", role=c("aut", "cre")), person("Thomas", "Lumley", role=c("ctb", "trl"), - comment="original S->R port and R maintainer until 2009")) + comment="original S->R port and R maintainer until 2009"), + person("Atkinson", "Elizabeth", role="ctb"), + person("Crowson", "Cynthia", role="ctb")) Description: Contains the core survival analysis routines, including definition of Surv objects, Kaplan-Meier and Aalen-Johansen (multi-state) curves, Cox models, @@ -21,9 +22,11 @@ License: LGPL (>= 2) URL: https://github.com/therneau/survival NeedsCompilation: yes -Packaged: 2019-12-01 17:42:21 UTC; therneau +Packaged: 2020-03-05 18:24:20 UTC; therneau Author: Terry M Therneau [aut, cre], Thomas Lumley [ctb, trl] (original S->R port and R maintainer until - 2009) + 2009), + Atkinson Elizabeth [ctb], + Crowson Cynthia [ctb] Repository: CRAN -Date/Publication: 2019-12-03 17:30:02 UTC +Date/Publication: 2020-03-07 12:10:03 UTC diff -Nru survival-3.1-8/inst/CITATION survival-3.1-11/inst/CITATION --- survival-3.1-8/inst/CITATION 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/inst/CITATION 2020-02-10 16:11:13.000000000 +0000 @@ -1,8 +1,11 @@ +year <- sub("-.*", "", meta$Date) +note <- sprintf("R package version %s", meta$Version) + bibentry(bibtype="Manual", - title = "A Package for Survival Analysis in S", + title = "A Package for Survival Analysis in R", author= person(c("Terry M"), "Therneau"), - year =2015, - note ="version 2.38", + year = year, + note = note, url="https://CRAN.R-project.org/package=survival", key= "survival-package" ) Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/adjcurve.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/adjcurve.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/approximate.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/approximate.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/compete.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/compete.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/concordance.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/concordance.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/multi.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/multi.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/other.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/other.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/population.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/population.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/splines.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/splines.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/survival.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/survival.pdf differ diff -Nru survival-3.1-8/inst/doc/survival.Rnw survival-3.1-11/inst/doc/survival.Rnw --- survival-3.1-8/inst/doc/survival.Rnw 2019-10-28 15:08:55.000000000 +0000 +++ survival-3.1-11/inst/doc/survival.Rnw 2020-01-28 21:37:52.000000000 +0000 @@ -460,10 +460,10 @@ Other options: \begin{itemize} - \item xaxt('r') It has been traditional to have survival curves touch + \item xaxs('r') It has been traditional to have survival curves touch the left axis (I will not speculate as to why). - This can be accomplished using \code{xaxt='S'}, which was the default in + This can be accomplished using \code{xaxs='S'}, which was the default before survival 3.x. The current default is the standard R style, which leaves space between the curve and the axis. \item The follow-up time in the data set is in days. This is very common in Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/tiedtimes.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/tiedtimes.pdf differ Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/timedep.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/timedep.pdf differ diff -Nru survival-3.1-8/inst/doc/timedep.R survival-3.1-11/inst/doc/timedep.R --- survival-3.1-8/inst/doc/timedep.R 2019-12-01 17:42:19.000000000 +0000 +++ survival-3.1-11/inst/doc/timedep.R 2020-03-05 18:24:18.000000000 +0000 @@ -125,20 +125,20 @@ ################################################### -### code chunk number 11: timedep.Rnw:583-584 +### code chunk number 11: timedep.Rnw:585-586 ################################################### attr(pbc2, "tcount") ################################################### -### code chunk number 12: timedep.Rnw:586-588 +### code chunk number 12: timedep.Rnw:588-590 ################################################### #grab a couple of numbers for the paragraph below atemp <- attr(pbc2, "tcount")[2:3,] ################################################### -### code chunk number 13: timedep.Rnw:669-675 (eval = FALSE) +### code chunk number 13: timedep.Rnw:671-677 (eval = FALSE) ################################################### ## temp <- subset(pbc, id <= 312, select=c(id:sex, stage)) ## pbc2 <- tmerge(temp, temp, id=id, death = event(time, status)) @@ -316,7 +316,7 @@ ################################################### -### code chunk number 28: timedep.Rnw:1176-1183 +### code chunk number 28: timedep.Rnw:1178-1185 ################################################### function(x, t, riskset, weights){ obrien <- function(x) { @@ -328,7 +328,7 @@ ################################################### -### code chunk number 29: timedep.Rnw:1193-1195 +### code chunk number 29: timedep.Rnw:1195-1197 ################################################### function(x, t, riskset, weights) unlist(tapply(x, riskset, rank)) diff -Nru survival-3.1-8/inst/doc/timedep.Rnw survival-3.1-11/inst/doc/timedep.Rnw --- survival-3.1-8/inst/doc/timedep.Rnw 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/inst/doc/timedep.Rnw 2020-02-06 17:04:24.000000000 +0000 @@ -278,7 +278,7 @@ taken from \code{data2}. The idea behind the function is that each addition will be ``slipped in'' to the original data in the same way that one -would slide a new card into an existing deck of cards. +would add a new folder into a file cabinet. It is a complex function, and we illustrate it below with a set of examples that sequentially reveal its features. @@ -326,7 +326,7 @@ \begin{itemize} \item Each row of the resulting data set represents a time interval (time1, time2] which is open on the left and closed on the right. - Covariate values for that row are the covariate values that apply + Covariate values for each row are the covariate values that apply over that interval. \item The event variable for each row $i$ is 1 if the time interval ends with an event and 0 otherwise. @@ -417,8 +417,8 @@ \subsection{Stanford heart transplant} The \code{jasa} data set contains information from the Stanford heart -transplant study, in the form that it appeared in the -paper of Crowley and Hu \cite{Crowley77}. +transplant study, in the form that it appeared in the J American Statistical +Association paper of Crowley and Hu \cite{Crowley77}. The data set has one line per subject which contains the baseline covariates along with dates of enrollment, transplant, and death or last follow-up. @@ -439,7 +439,7 @@ treatment for days 1--10 and as transplanted starting on day 11. That is, except for patient 38 who died on the same day as their procedure. They should be treated as a transplant death; the problem is resolved by - moving this transplant back .5 day. + moving this forward in time by .5 day. \item The treatment coefficients in table 6.1 of the definitive analysis found in @@ -497,7 +497,9 @@ \code{tdc} or \code{cumtdc} call that has two arguments and 0 in all other cases. If the variable being created is already a part of data1, -then our updates make changes to that variable. +then our updates make changes to that variable; in the case of +an \code{event} call it results in an addition, for \code{tdc} +it results in wholesale replacement of any prior values. Be careful of this. This feature is what allowed for the \code{infection} indicator to be build up incrementally in the cgd example above, but quite surprising results can occur when you Binary files /tmp/tmpmVqWwN/ImwY0RRspX/survival-3.1-8/inst/doc/validate.pdf and /tmp/tmpmVqWwN/W1H07LbZYR/survival-3.1-11/inst/doc/validate.pdf differ diff -Nru survival-3.1-8/inst/NEWS.Rd survival-3.1-11/inst/NEWS.Rd --- survival-3.1-8/inst/NEWS.Rd 2019-12-01 16:48:18.000000000 +0000 +++ survival-3.1-11/inst/NEWS.Rd 2020-03-05 13:29:54.000000000 +0000 @@ -1,8 +1,38 @@ \name{NEWS} \title{NEWS file for the survival package} +\section{Changes in version 3.1-11}{ + \itemize{ + \item Add 'n' vector to survcheck, update tests. + \item Improve multistate coxph printout +}} + +\section{Changes in version 3.1-10}{ + \itemize{ + \item Fix error in concordance(fit, newdata= xxx): it used the old + data at one point. Pointed out by Matjin van't Hooft. + \item Minor update to tmerge -- print a warning when a user attempts + to update a variable using tdc, namely that this results in complete + replacement of said variable, with no retention of prior values. + \item Minor fixes to plot.survfit, when both xlim and xaxs='S' are + specified. Error and suggested fixes from Heinz Tuechler. + \item Two failures in makepredictcall.pspline, found by 2 different + users, both cases of options that should have been carried forward. + A subsequent termplot() call failed if pspline had 'nterm' or + 'combine'. + \item Fix memory leak (forgotten PROTECT) in zph[12].c; caught by B Ripley +}} + +\section{Changes in version 3.1-9}{ + \itemize{ + \item The plot.cox.zph function would fail for a stratified model, + if some covariate was a constant in one of the strata. The cox.zph + routine correctly labeled those points as uninformative via NA, the + plot routine needed to actually ignore them. +}} + \section{Changes in version 3.1-8}{ \itemize{ - \item For a multistate model which no covariates, the start.time + \item For a multistate model with no covariates, the start.time value (if present) was inadvertently omitted from the result. This affects the default plot. \item Per user request, the tmerge function can now use Dates as the diff -Nru survival-3.1-8/man/coxph.Rd survival-3.1-11/man/coxph.Rd --- survival-3.1-8/man/coxph.Rd 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/man/coxph.Rd 2020-02-25 18:20:26.000000000 +0000 @@ -91,8 +91,8 @@ the purposes of a robust variance. If present, it implies \code{robust}. This variable will normally be found in \code{data}.} - \item{istate}{optional variable giving the intial state for each - subject. This variable will normally be found in \code{data}.} + \item{istate}{optional variable giving the current state at the start + each interval. This variable will normally be found in \code{data}.} \item{statedata}{optional data set used to describe multistate models.} \item{model}{ logical value: if \code{TRUE}, the model frame is returned in component @@ -185,7 +185,24 @@ this case the \code{tt} argument will be a function or a list of functions (if there are more than one tt() term in the model) giving the appropriate transform. See the examples below. + +One user mistake that has recently arisen is to slavishly follow the +advice of some coding guides and prepend \code{survival::} onto +everthing, including the special terms, e.g., + \code{survival::coxph(survival:Surv(time, status) ~ age + + survival::cluster(inst), data=lung)} +First, this is unnecessary: arguments within the \code{coxph} call will +be evaluated within the survival namespace, so another package's Surv or +cluster function would not be noticed. +(Full qualification of the coxph call itself may be protective, however.) +Second, and more importantly, the call just above will not give the +correct answer. The specials are recognized by their name, and +`survival::cluster' is not the same as `cluster'; the above model would +treat \code{inst} as an ordinary variable. +A similar issue arises from using \code{stats::offset} as a term, in +either survival or glm models. } + \section{Convergence}{ In certain data cases the actual MLE estimate of a coefficient is infinity, e.g., a dichotomous variable where one of the @@ -195,10 +212,10 @@ becomes effectively singular, an argument to exp becomes too large for the computer hardware, or the maximum number of interactions is exceeded. -(Nearly always the first occurs.) +(Most often number 1 is the first to occur.) The routine attempts to detect when this has happened, not always successfully. -The primary consequence for he user is that the Wald statistic = +The primary consequence for the user is that the Wald statistic = coefficient/se(coefficient) is not valid in this case and should be ignored; the likelihood ratio and score tests remain valid however. } @@ -228,7 +245,7 @@ long. With (start, stop) data it is much worse since the recursion needs to start anew for each unique start time. -A second issue is that of artificial ties due to floating-point +A separate issue is that of artificial ties due to floating-point imprecision. See the vignette on this topic for a full explanation or the \code{timefix} option in \code{coxph.control}. Users may need to add \code{timefix=FALSE} for simulated data sets. diff -Nru survival-3.1-8/man/survcheck.Rd survival-3.1-11/man/survcheck.Rd --- survival-3.1-8/man/survcheck.Rd 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/man/survcheck.Rd 2020-02-25 17:16:25.000000000 +0000 @@ -13,10 +13,10 @@ \item{data}{data frame in which to find the \code{id}, \code{istate} and formula variables} \item{id}{an identifier that labels unique subjects} - \item{istate}{an optional vector giving the initial state for each - observation, i.e., the state at the start of each interval} - \item{istate0}{default label for the initial state when \code{istate} - is missing} + \item{istate}{an optional vector giving the current state at the start + of each interval} + \item{istate0}{default label for the initial state of each subject (at + their first interval) when \code{istate} is missing} \item{timefix}{process times through the \code{aeqSurv} function to eliminate potential roundoff issues.} \item{\ldots}{other arguments, which are ignored (but won't give an @@ -56,7 +56,8 @@ \item{statecount}{table of the number of visits per state, e.g., 18 subjects had 2 visits to the "infection" state} \item{flags}{a vector giving the counts of each check} - \item{istate}{a revised initial state vector} + \item{istate}{a copy of the istate vector, if it was supplied; + otherwise a constructed istate that satisfies all the checks} \item{overlap}{a list with the row number and id of overlaps (not present if there are no overlaps)} \item{gaps}{a list with the row number and id of gaps (not present if @@ -66,7 +67,26 @@ \item{jumps}{a list with the row number and id of jumps (not present if there are no jumps} } -\author{Terry Therneau and Beth Atkinson} +\note{ +For data sets with time-dependent covariates, a given subject will often +have intermediate rows with a status of `no event at this time'. (numeric +value of 0). +For instance a subject who started in state 1 at time 0, transitioned to state +2 at time 10, had a covariate \code{x} change from 135 to 156 at time +20, and a final transition to state 3 at time 30. +The response would be \code{Surv(c(0, 10, 20), c(10, 20, 30), c(2,0,3))}: +the status variable records \emph{changes} in state, and there was no +change at time 20. +The \code{istate} variable would be (1, 2, 2); it contains the \emph{current} +state, and so the value is unchanged when status = censored. +Thus, when there are intermediate observations \code{istate} is not +simply a lagged version of the status, and may be more challenging to +create. +One approach is to let \code{survcheck} do the work: call it with +an \code{istate} argument that is correct for the first row of each +subject, or no \code{istate} argument at all, and then insert the +returned value into ones data frame. +} \keyword{ survival } diff -Nru survival-3.1-8/MD5 survival-3.1-11/MD5 --- survival-3.1-8/MD5 2019-12-03 17:30:02.000000000 +0000 +++ survival-3.1-11/MD5 2020-03-07 12:10:03.000000000 +0000 @@ -1,4 +1,4 @@ -71ec2bdee8ec27be0ac3d5ffcf7085d9 *DESCRIPTION +3b697367d7d7d58342630d0e9c27761e *DESCRIPTION 0821b94712becd7c2a2605a0dde65a04 *NAMESPACE 8832a09384db2ba7d905016e7b1e89ca *R/Surv.R 18d616dcd94d2763aa6d98a396f26361 *R/aareg.R @@ -18,12 +18,12 @@ 4696548963d14529a80b67935ec99dae *R/cipoisson.R 6a7ae6876a34e3f55d34de654dc3afe6 *R/clogit.R 142abff209b16e6fb74a2962442ffc34 *R/cluster.R -31bd032a14dff02017b2d6df4c8949cf *R/concordance.R +c7a7706c0c573af3532ddf7055ea1663 *R/concordance.R 18f979008cc201d4d1637eadb704076e *R/cox.zph.R c35eeb88e9603b732863591821b83cf9 *R/coxexact.fit.R 94e49b45d947935252dac8935c85b14f *R/coxpenal.df.R 17d8987218ae14a72cf60a07338df2e1 *R/coxpenal.fit.R -daba44cb14d4dd74a4dc3b431e356886 *R/coxph.R +1e413d748658755fc41396f8121775c4 *R/coxph.R 1db9683e308bfaba24c8af91d9a509d7 *R/coxph.control.R fbe682c59326ff0b09748658e83a7477 *R/coxph.detail.R 62c4314464ee0a2a66e880464c913997 *R/coxph.fit.R @@ -59,31 +59,31 @@ 2e8efcaf679a4cabbdaec2b342d38d51 *R/pade.R 00367286e34dd8b6e45ec7b326c0e355 *R/parsecovar.R 5220681f176b87fb822910dea854ccf4 *R/plot.aareg.R -bfbeb021ec253a85a848e0c863c2bae4 *R/plot.cox.zph.R -7cc823f60a9521b6d729cb33b17e1060 *R/plot.survfit.R +dfc82fb9e8f44bbe7cafe8ca9a02f09e *R/plot.cox.zph.R +e9cd2997bb6db5e524e0a9770f9c3676 *R/plot.survfit.R 7a866d546f5264744430536c779d76ed *R/predict.coxph.R bf8f5ab0bb5491ea6ab6c1c3ac4fdb5f *R/predict.coxph.penal.R ab222d970b26ee908f2acc250b04e0fe *R/predict.survreg.R a2e60c0b4dd9667bcfbac613c929b42c *R/predict.survreg.penal.R e8dc0601452488a22df48d545be06edc *R/print.aareg.R -98adb33cab7deccf635903ab243eaf8d *R/print.coxph.R +99ebd8fbe5183ec676ba7998ca82f349 *R/print.coxph.R 9bc844ff207051b2c422f17cfb026981 *R/print.coxph.null.R f1c8f318572434c69c5ec8fff7b6c108 *R/print.coxph.penal.R 4b59fc74eb3ce7c6dd195651339a09d6 *R/print.pyears.R 4e794cd106b6225d687ed58acfe3c112 *R/print.ratetable.R -ab982a110a0bf50fef3dee4d42ebc5b2 *R/print.summary.coxph.R +b8493f405943a8c2f3191af9033dc894 *R/print.summary.coxph.R 592a17659e17a6ccb85d519c4b284034 *R/print.summary.coxph.penal.R a5fbce4916ba3ac805bf29e434eedf06 *R/print.summary.survexp.R fa8753120300851242ec7286d3edfa64 *R/print.summary.survfit.R 0823b87d28b7b14ae423d8b92c1c63ad *R/print.summary.survfitms.R 69cfc284bb126050a4adbb0f2aeddd50 *R/print.summary.survreg.R -f34d207de2825bd5139f4a76c194f385 *R/print.survcheck.R +9f0a3fa3712d4b8beb8fc2822e341915 *R/print.survcheck.R 69c1cd6a55cca3dfad2e1d65745a5692 *R/print.survdiff.R 435591842ba7380ea9c94033fd867165 *R/print.survexp.R 90f112bfdd842ec827f41b268381e0e7 *R/print.survfit.R 67edc10ae58de66955889a85822297b7 *R/print.survreg.R 4227591baa8a237ff83b8976b230856b *R/print.survreg.penal.R -0a6b806244056b84905449d75758979d *R/pspline.R +fad4a4413fb7af9cbf07a3445d3bd3e0 *R/pspline.R dd5c3a90d9adced8b88334d5fbed68cc *R/pyears.R c596219c1d6f31581a3ec622c586b0f7 *R/quantile.survfit.R b2f82c369fd6c8bf4496b1fe56ddef5d *R/ratetable.R @@ -92,14 +92,14 @@ cc274d54216721d3964568f28b2f1e21 *R/residuals.coxph.R 871dc4bdd8fc9c0528e9c794498afcab *R/residuals.coxph.null.R 71bb8b32a975d5add7bff7d5963c8ea2 *R/residuals.coxph.penal.R -754627902dd32d5dc894506a03aedbcb *R/residuals.survreg.R +3bd54c45e897334ffd3360e1c5a447f8 *R/residuals.survreg.R 754a417ad23f70056005558688c57e88 *R/residuals.survreg.penal.R 5a28f06bc0463f5278e17e894b5bbc78 *R/ridge.R 1708f9c9e6bc1d7fcb9f209f47f5fe5e *R/stacker.R c1bc3f2ca2758374d291aae108e9448a *R/statefig.R d4965488f4b863ef65b00c93886454e9 *R/strata.R cf9ec91bb29089a64cae0e83801e0b90 *R/summary.aareg.R -fa194b5d7d4908607fcb2321b5a4d033 *R/summary.coxph.R +0c2e3984f9abfe7a25c54cf90f418ca3 *R/summary.coxph.R ab7c37261a9ad2ea16b24155c4b75c41 *R/summary.coxph.penal.R 50a7f4971ccb0004da11691bd7dbd9c9 *R/summary.ratetable.R 5c340d79d336b2033f2287b104bfffd8 *R/summary.survexp.R @@ -108,7 +108,7 @@ 3b53d1f647391695c4bac97f0ab201b8 *R/survConcordance.fit.R d77d5fb984e1cd3345e55df7af866ef9 *R/survSplit.R bc393f42ff473330532bbdbf046ea772 *R/survcallback.R -5c48b721ba5beef18554b6df5c63f358 *R/survcheck.R +b9a030c232021e484ee1d316e0bc3851 *R/survcheck.R 82c4bc1c9eb611f6350036873b32431b *R/survdiff.R c6ddf6a691e8b761df6e38e5bb6f8884 *R/survdiff.fit.R b045d10cf20e99cab84090c2efcd1d6a *R/survexp.R @@ -132,7 +132,7 @@ 7e396c0fd8c8b7b97f5e5dcf79064396 *R/survreg.fit.R 7c797d7ee0fcfd0a92dd654a7e5aa03d *R/survregDtest.R 016e8b5dd0f4dfe6fb13b4f262bec247 *R/tcut.R -6636841371b27dc4733bf0803f80a1eb *R/tmerge.R +f18e4b81f3bcc1861ee866609448656c *R/tmerge.R 5bd109e83e16ef7659f0a81709694ec2 *R/untangle.specials.R 6a0998da3982b925dcb70fa95978dc36 *R/xtras.R 0e76b8a5afc0c8207059775c65423783 *R/yates.R @@ -166,43 +166,43 @@ 6b016a6be8a1ff47f01e307f8e5734a8 *data/udca.rda 76fb141b816ea857308105f26e26e87f *data/uspop2.rda d4b78dd46512d982683f6bb2a025aa9a *data/veteran.rda -b6991f4124dcc3bef51e7340d91beb64 *inst/CITATION +ba01ac7a9482cc481afa2bcfd8126e99 *inst/CITATION 3e8ebae349953a3b7a665eecce562802 *inst/COPYRIGHTS -fbb0fac5984cb26dc938cc24f39530fa *inst/NEWS.Rd +e924d8b407ebf63bf6dd3ac27fdcbbc7 *inst/NEWS.Rd 097160093b374a96be83fb1c1f6265c9 *inst/doc/adjcurve.R eb3426c1018a1a407180b922c8568920 *inst/doc/adjcurve.Rnw -db65ddbe69c6c4f1c07238d7f9c3690e *inst/doc/adjcurve.pdf +e54b0052d616b7369ca6eca2f36d82d3 *inst/doc/adjcurve.pdf ae5a9072b9763462ae72b58e0b15d38f *inst/doc/approximate.R 0f1ca788773dcad8ef5f1372722e629c *inst/doc/approximate.Rnw -7451c1da808b858024ad872b8814fb4c *inst/doc/approximate.pdf +c25a20627f776b086f3a43df54d775a0 *inst/doc/approximate.pdf f5db3b691ac914fbf44a096ec5b43470 *inst/doc/compete.R 2efc66310e74010eb99efcc6cc4c6fa1 *inst/doc/compete.Rnw -ee3e6d28c772581414892a5067d3c73e *inst/doc/compete.pdf +ea4edd5522d10ea5c8a1952c1334819c *inst/doc/compete.pdf d263e16c5004ad8e0b718d8deee81a29 *inst/doc/concordance.R b931633a7cacd56bde2ecf4654cebc96 *inst/doc/concordance.Rnw -0b82d3a871b2bd9e74434c599b4c8a21 *inst/doc/concordance.pdf +03d7eeef6380ca909c7a35458ce04b4a *inst/doc/concordance.pdf 5f0a233bbc83852d406f28858d6cacdb *inst/doc/multi.Rnw -80829a6d825dbd2b70d11abb327967ce *inst/doc/multi.pdf +d33fb973fa8846c6cb2022b9cf451759 *inst/doc/multi.pdf 24c8cf5cb39d43041713e77c4d88d6a0 *inst/doc/other.Rnw -66392c947560913a07441242b8c1c00b *inst/doc/other.pdf +ae36e03a92fb0b2ab487cd9835904044 *inst/doc/other.pdf 24c8009f8a654c9d1bde3a8371385a11 *inst/doc/population.R e5e28360d2f5718b176951da04d5e24b *inst/doc/population.Rnw -ad74988aa5265e5d5908254bdeb7d1b9 *inst/doc/population.pdf +a336982b669cb4a5ce50d2d36096a7ee *inst/doc/population.pdf 3b2499552d20bd20231bca68696caf23 *inst/doc/splines.R d05b97e79328012602fd6e7a2867ed81 *inst/doc/splines.Rnw -0eef2363060500a83f0654f8f9b581d2 *inst/doc/splines.pdf +98589c1b235d5b6e41d31707fc1eff9a *inst/doc/splines.pdf da96690e64daf512296d6d99dcbb25a6 *inst/doc/survival.R -373c53ef508d72a258a552759e7d42fe *inst/doc/survival.Rnw -eb4ffdc4d2c51524c0d199b8bd0618a9 *inst/doc/survival.pdf +1d7f54a86d33cea87b3724403d4771e5 *inst/doc/survival.Rnw +8455c77e9073fadc0763ef98bd97be75 *inst/doc/survival.pdf 40ceb0c930e083a97a0a35702af6c5da *inst/doc/tiedtimes.R e786486fd295208ebdc6a15d3fe56e5b *inst/doc/tiedtimes.Rnw -1fbfb8d5ead888bf36923d77b8598234 *inst/doc/tiedtimes.pdf -c9caef55c6e902140899b994c35ed1a6 *inst/doc/timedep.R -5d034a0848ea5caa82a65639264b9782 *inst/doc/timedep.Rnw -e69c583da6383d06e989bd838d407a34 *inst/doc/timedep.pdf +00127212c1d88d917cbb83fe13fe1b9f *inst/doc/tiedtimes.pdf +1f618f32026fae363a561283857dc412 *inst/doc/timedep.R +a26837c352477a1d288bad00252022d9 *inst/doc/timedep.Rnw +b0d6a589d1f886ee893cbc09a8d15e58 *inst/doc/timedep.pdf 30e95a3026ae131c606f45fdb4732a4f *inst/doc/validate.R dfffceef39872c105be8fb3b54d2c15d *inst/doc/validate.Rnw -1b24a81d4b20debf69ddb797a5a4f784 *inst/doc/validate.pdf +56210604d24484e3e48055e11dd43a59 *inst/doc/validate.pdf 35602740015d6a49b48756564c05c191 *man/Surv.Rd bffb4bc6a2827e87f8a1b38d0cf3c083 *man/Survmethods.Rd 2712884a7d58a481678f80f0ff17a0a9 *man/aareg.Rd @@ -223,7 +223,7 @@ ad5e45c483b88a6f5476e26427538192 *man/concordance.Rd b3d81c8e8d7219c0588845a5019cb213 *man/concordancefit.Rd 1f91aedb37f89d1951ded0d47a4a0c7b *man/cox.zph.Rd -faeff07a36cd0bafbd858617f7c5ea64 *man/coxph.Rd +e4ac334cb8f28e44bd70ca0c9e183d4b *man/coxph.Rd 3aa398cfa6989e15a73eef4cd65fe565 *man/coxph.control.Rd d90b227d1ff0e2c44c2382944f05ddb6 *man/coxph.detail.Rd ba16ad811717398eb585ea042731939e *man/coxph.object.Rd @@ -288,7 +288,7 @@ 4df84fa21a9d207959343b48b0009ee8 *man/summary.survfit.Rd 84ee7e8d8a38c7651b2882373dcb0f32 *man/survConcordance-deprecated.Rd 54506d288b5871c1d8588e3e80df3ce9 *man/survSplit.Rd -5417f895c80be70c93fde7853d993b75 *man/survcheck.Rd +68250ff6c2035fc4a3234f0eff598b35 *man/survcheck.Rd d3e26ce68b91b1ea2facbcf3d5197b5b *man/survdiff.Rd f1d6904a6d2c6c38493f9b928e15ad06 *man/survexp.Rd 65e8c2fdd454ccb1a16f6dcdb86d77bc *man/survexp.fit.Rd @@ -325,22 +325,23 @@ 87b794e23fcfbe3c21e73615ea7cbb97 *noweb/agreg.Rnw 064c19c6153c5934c977c3e06ed84df9 *noweb/casecohort.Rnw deee274f667f087bb948bdee592d076b *noweb/cch.Rnw -5c0e1f2f964e48d94e428611d0c2b7fb *noweb/code.nw -8610f3b2df8e5559b151317a8119082f *noweb/code.toc -3a3d3aee3a4a9a7b19a5f0c45fd6551f *noweb/concordance.Rnw +a43db188eefc16713134fe268b5e9bda *noweb/code.nw +3f905b0f2df03a5142ef3433f9ba05f1 *noweb/code.toc +b735fea34a1890db2a750a5bd029599c *noweb/concordance.Rnw 099ea7b709d6521fc862938c279551f5 *noweb/concordance2.Rnw -0c7d9cd6b8ef7322606fba80f3e35999 *noweb/coxph.Rnw -0b03327c96b83acec21eb57563b151b4 *noweb/coxph.Rnw.orig +e1fa07b0744e1d47e3e3e45965c523c9 *noweb/coxph.Rnw +57afceed6073aebc91bb02f76c32217e *noweb/coxph.Rnw.orig 97e9e74a285071aaac11f8c36ae754f7 *noweb/coxsurv.Rnw e2d0b6ad8a1886a63b56d0362a8eed2d *noweb/coxsurv2.Rnw 8a445a83e715f284ff318d1c6bde722c *noweb/coxsurv3.Rnw 01b1d864a903a3c3cff1ea0725075363 *noweb/exact.nw dd53136ac036954c3ac43c1f7d95050f *noweb/finegray.Rnw ec7ddcb0dd1cb1bfd3a3722d2da12cfc *noweb/main.Rnw -31aad4320b375a1c097ec48087e4589f *noweb/msurv.nw +696fee2a1f807fe3f01ef5d4f86db186 *noweb/maximacode +a515e93036e8a44d9d1a5b4ed8b4d9d9 *noweb/msurv.nw df28925c86bc02dbd9bf8b2f73856f4f *noweb/noweb.sty 1284d9fb78ba492253ccceefea18b17b *noweb/parse.Rnw -fbafe763beacc125abdc6fc11cf0892b *noweb/plot.Rnw +cd8fbfc78d8c4ba12e8203c38800955c *noweb/plot.Rnw 354522d4762cf4cac3189ae1404ce6cf *noweb/predict.coxph.Rnw 4c349d0baeb6aa069e313feb6ff01fc9 *noweb/pyears.Rnw 36cf079089d42e45b7a7f4ef6ea7755f *noweb/pyears2.Rnw @@ -362,7 +363,7 @@ 8337d1debe3db49ad09a3ab8c9872789 *noweb/rates/usinfant.dat 4cea128faf7fbd8e6ab46070dce67cb3 *noweb/ratetable.Rnw 00153a591d69dc50d6852fb8f9f2ddb2 *noweb/refer.bib -fd02d5e7a239e8ebcb06296bbdd22167 *noweb/residuals.survreg.Rnw +7609336443d38f08fd88c62f70bf3953 *noweb/residuals.survreg.Rnw c7315b78224994c56c4d140e221fab3a *noweb/statefig.Rnw 67c10c8c39f68add5d487185eb3f9f7e *noweb/survConcordance.Rnw 2e1d48d5dbfdc85b69dae483686173e4 *noweb/survexp.Rnw @@ -371,7 +372,8 @@ 64198344a2855e8c6d3f78e09a3968c6 *noweb/survfitKM.Rnw 9e4f3d256bf768545064ac2b8044d361 *noweb/survfitms.Rnw 9534b235ab22d3c6d25f7a1821cac049 *noweb/tail -2d864277a77464462f4d5360d7334343 *noweb/tmerge.Rnw +cd92b1bc6874cf319aaf31582830e052 *noweb/tmerge.Rnw +6db146f5a9dd7b8760d81c0416d86897 *noweb/tmerge.Rnw.orig ba373139f6ca0bee540a95f16f65e9e3 *noweb/yates.Rnw 7b4c8ccf7679d96d4a13e659b08116f8 *noweb/yates2.Rnw 9ae7d3e5048dd95d39a5a7210bad5eb0 *noweb/zph.Rnw @@ -412,8 +414,8 @@ 6a5de3098ea717defac73b46aaf36006 *src/doloop.c 893b8f26f882396c490b9c18b3072051 *src/finegray.c d581082f66939c0281ce8c859f723148 *src/gchol.c -5baa10a3c202f9ada646b0380822db30 *src/init.c -9d4c35fa01498378f366900acce90441 *src/multicheck.c +779982ba171c2f19f696a2ca06aeb18d *src/init.c +804f9cb3e2b49a270b6b011ea14f5ac9 *src/multicheck.c 1385cc1bf5e1b1de23bd1a842be745d4 *src/norisk.c 23bb38ccf2c88c0f47776f5b38b03ab6 *src/pyears1.c 3f85ef9d57e28eb78ac14afb2cf8e855 *src/pyears2.c @@ -426,15 +428,15 @@ 8dc673be20ef1e7ad592d053336a5ea3 *src/survfitci.c a58fdee2ebab7f4c5e13b982c72868ad *src/survfitkm.c 0ca9cc1702f9f87b2e25c414892f4e14 *src/survpenal.c -75a12a1367acdd8058a910b9168f8d73 *src/survproto.h +4de707697ae1a81cfa1b1f4b5fbeb678 *src/survproto.h 2f1657fe25d2d8d6c1ecd3ae64b8a658 *src/survreg6.c 7483bb0e3da8485fae01944322021385 *src/survreg7.c fec3ad6f2d4918cb7115f6efe5439bf3 *src/survregc1.c 803c4db11e03621954ac362de4957b48 *src/survregc2.c 1f287cced7b7527881c56b5e3f7c5a70 *src/survsplit.c fd828f24cd4f2e4b68bd77699c5b68a6 *src/tmerge.c -a5902a7b6c356734387bf0e17d349904 *src/zph1.c -e4a7135b73d591c55444f6265cf761c4 *src/zph2.c +74f6ac54c5dc76cda6575cd88724e98f *src/zph1.c +a86b40542fa5e2207d7b527e67aeba6c *src/zph2.c ab2989d395c2102a8802f0734419eb22 *tests/aareg.R b9c7ea3cc2e4ac2607c9c79e495b1d51 *tests/aareg.Rout.save 7adcf3da61f0cb888ed3ae1787d85fd5 *tests/anova.R @@ -450,9 +452,9 @@ e2dac82a95f6c623bc4685b1b9c4c389 *tests/book4.R 6187357af28e99b9a2ecf99b06385f86 *tests/book4.Rout.save b878fcc994199cbc85e659c193266496 *tests/book5.R -bf8c822bd6fc2b2a47a9892b26d9193d *tests/book5.Rout.save +60a26d189bc9c0a52cd44995289638d6 *tests/book5.Rout.save 6c123e7a51c707a4708de781550425a8 *tests/book6.R -08d7daf71984ff2501e0be471c51921f *tests/book6.Rout.save +bebf1ddbc1915a05435632589db34e64 *tests/book6.Rout.save 7f0043933913df6795c66c8af6a278fc *tests/book7.R b3297d78ca2a2d942a91ab83954c8123 *tests/book7.Rout.save 047709f7cab960044afaefd038e09b1e *tests/cancer.R @@ -489,7 +491,7 @@ 2d9aa1e5509fbc7504b3f7d729e48911 *tests/doaml.R 6c7a563e04818392a55bfb6aec2e55db *tests/doaml.Rout.save cdf564fec0e0177ca97abbd83c70499a *tests/doweight.R -30e61555f74f20162832dda71a93e634 *tests/doweight.Rout.save +e7d0b517b5468eece32273c0526a702a *tests/doweight.Rout.save 106d5e11a7b70ae5712b1a631236b303 *tests/expected.R 067cf68e6507773929a7f0a76f405a45 *tests/expected.Rout.save d0e4b970ff89cca427af7fa83819529e *tests/expected2.R @@ -555,7 +557,7 @@ 672df5418ed229dc3e76949585e3790e *tests/r_resid.R 5687da81357331f70ecbde12b940b6a6 *tests/r_resid.Rout.save fa316196fb935eacc2d179388d728ce9 *tests/r_sas.R -1a3e481043d7aee83a08323da5435c9a *tests/r_sas.Rout.save +a0b832e16f783c84ee3e706a38c1c7a0 *tests/r_sas.Rout.save 59912fe4c23adf729630b2a9256749bd *tests/r_scale.R 9c44d08c1462edc432d1162f6353a785 *tests/r_scale.Rout.save 3c6a7056af6fb278f774abf6a69c3968 *tests/r_stanford.R @@ -568,6 +570,7 @@ bcf83472af25c444668fb004755e2b52 *tests/r_user.Rout.save 1f42efa1bfc2374f79ba2bc55ef9cfdb *tests/ratetable.R f33d097be853fbabb5536dc1332c5e90 *tests/ratetable.Rout.save +27d16e537f3d56a0d9a25b49830dc8e2 *tests/sexpm.save 5a65fb2b39182a19953db1807ec3f991 *tests/singtest.R 3a893b1ed2457350aed62c2cd17ab441 *tests/singtest.Rout.save f7a452292c0f25593b8cbc15bdd572d5 *tests/strata2.R @@ -580,8 +583,8 @@ ae2c73dff8c70a7e30d28ed7a5a8bee4 *tests/surv.Rout.save 85c94bf1652e4b3102d0841de108cd9e *tests/survSplit.R 38f51a5663deb639f278e410c1a26568 *tests/survSplit.Rout.save -9d699784b63cb37da438af411d85b50e *tests/survcheck.R -6078de8f6f7a35c9b7fc61db3bad0f58 *tests/survcheck.Rout.save +a64057b3067b7f26825fcb7c64061618 *tests/survcheck.R +22eb181aa50b34319531e039c8e4b314 *tests/survcheck.Rout.save 7740b37978d2b72f2f554f4cf42022f9 *tests/survfit1.R 1a0bc7dfb07c7d9d2092cb01a5c5e73e *tests/survfit1.Rout.save 41e3629b6abc984689331b047ee491ba *tests/survfit2.R @@ -607,8 +610,8 @@ ecefc6b3b922ee53c8d7234b519791d5 *tests/tmerge.Rout.save 6dbf149f8da192b5dbbbcf0a2cc47840 *tests/tmerge2.R 81ce8cb239643765af84eab1fb45f170 *tests/tmerge2.Rout.save -b293576661a786d575b375ca16d78ac9 *tests/tmerge3.R -eb2fd4977d5ac90e9f38804d3a53293d *tests/tmerge3.Rout.save +ba7b9439c287664fb97b0553ce59d06d *tests/tmerge3.R +fab7f4afd71506b6f2fb5efa3c328490 *tests/tmerge3.Rout.save 1e369715fdeb935e8cc784f0f40502ce *tests/tt.R 1bdb4feff5648ee13239625ef5d090f9 *tests/tt.Rout.save b1dd6e9b7bf7a5a26e8879c3f7fb03d9 *tests/turnbull.R @@ -631,7 +634,7 @@ e5e28360d2f5718b176951da04d5e24b *vignettes/population.Rnw 7d2afaa75cefcc21619b18780d9bb86b *vignettes/refer.bib d05b97e79328012602fd6e7a2867ed81 *vignettes/splines.Rnw -373c53ef508d72a258a552759e7d42fe *vignettes/survival.Rnw +1d7f54a86d33cea87b3724403d4771e5 *vignettes/survival.Rnw e786486fd295208ebdc6a15d3fe56e5b *vignettes/tiedtimes.Rnw -5d034a0848ea5caa82a65639264b9782 *vignettes/timedep.Rnw +a26837c352477a1d288bad00252022d9 *vignettes/timedep.Rnw dfffceef39872c105be8fb3b54d2c15d *vignettes/validate.Rnw diff -Nru survival-3.1-8/noweb/code.nw survival-3.1-11/noweb/code.nw --- survival-3.1-8/noweb/code.nw 2019-12-01 17:40:18.000000000 +0000 +++ survival-3.1-11/noweb/code.nw 2020-03-05 13:32:32.000000000 +0000 @@ -388,7 +388,7 @@ Those rows can be removed from the model frame before creating the X matrix. Also identify partially used rows, ones where the necessary covariates are present for some of the possible transitions but not all. -Those obs are dealt with later by the stackup function. +Those obs are dealt with later by the stacker function. <>= # first vector will be true if there is at least 1 transition for which all # covariates are present, second if there is at least 1 for which some are not @@ -435,15 +435,22 @@ Also create the initial values vector. The stacker function will create a separate block of observations for every -unique value in \code{cmat[1,]}. +unique value in \code{cmap[1,]}. Now say that two transitions A:B and A:C share the same baseline hazard. Then either a B or a C outcome will be an ``event'' in that stratum; they would only be distinguished by perhaps having different covariates. +The first thing we do with the result is to rebuild the transitions matrix: +the current version was created before removing missings, and can +seriously overstate the number of transitions available. +Then set up the data. <>= cmap <- parsecovar3(tmap, colnames(X), attr(X, "assign")) xstack <- stacker(cmap, as.integer(istate), X, Y, strata=istrat, states=states) +rkeep <- unique(xstack$rindex) +transitions <- survcheck2(Y[rkeep,], id[rkeep], istate[rkeep])$transitions + X <- xstack$X Y <- xstack$Y istrat <- xstack$strata @@ -5500,7 +5507,7 @@ else stop("left hand side of the formula must be a numeric vector or a survival object") } if (timefix) y <- aeqSurv(y) - rval <- list(y= y, x= predict(object)) + rval <- list(y= y, x= predict(object, newdata)) # the type of prediction does not matter, as long as it is a # monotone transform of the linear predictor } @@ -5768,7 +5775,7 @@ if (is.null(wt)) vmat <- crossprod(dfbeta) else vmat <- t(wt * dfbeta) %*% dfbeta rval <- list(concordance=concordance, count=count, - n=length(flist[[1]]$x), var=vmat, + n=flist[[1]]$n, var=vmat, cvar= sapply(flist, function(x) x$cvar)) if (influence==1) rval$dfbeta <- dfbeta @@ -5917,7 +5924,7 @@ Referring again to figure \ref{treefig}, \code{btree(13)} yields the vector \code{7 3 8 1 9 4 10 0 11 5 12 2 6} meaning that the smallest element -will be in position 8 of the tree, the next smallest in position 2, etc, +will be in position 8 of the tree, the next smallest in position 4, etc, and using indexing that starts at 0 since the results will be passed to a C routine. The code just above takes care to do all arithmetic as integer. @@ -8619,7 +8626,9 @@ else score <- cbind(score, deriv[,4]) } rr <- score %*% vv - if (type=='dfbetas') rr <- rr %*% diag(1/sqrt(diag(vv))) + # cause column names to be retained + # old: if (type=='dfbetas') rr[] <- rr %*% diag(1/sqrt(diag(vv))) + if (type=='dfbetas') rr <- rr * rep(1/sqrt(diag(vv)), each=nrow(rr)) if (type=='ldcase') rr<- rowSums(rr*score) } @@ -10751,7 +10760,7 @@ @ \subsubsection{C-code} (This is set up as a separate file in the source code directory since -it is easier to make the code stay in C-mode if the file has a .nw +it is easier to make emacs stay in C-mode if the file has a .nw extension.) <>= @@ -12070,7 +12079,7 @@ else stop("invalid cumhaz argument") } else if (inherits(x, "survfitms")) { - i <- (x$states != noplot) + i <- !(x$states %in% noplot) if (all(i) || !any(i)) { # the !any is a failsafe, in case none are kept then ignore noplot ssurv <- smat(x$pstate) @@ -12256,7 +12265,7 @@ if (!missing(xlim) && !is.null(xlim)) { tempx <- xlim - if (xaxs == 'S') tempx[2] <- tempx[2] + diff(tempx)*1.04 + if (xaxs == 'S') tempx[2] <- tempx[1] + diff(tempx)*1.04 } else { temp <- stime[is.finite(stime)] @@ -12327,7 +12336,7 @@ if (resetclip) { # yes, do it - if (xaxs=='S') tempx <- c(tempx[1], temp[1]) + #if (xaxs=='S') tempx <- c(tempx[1], temp[1]) clip(tempx[1], tempx[2], tempy[1], tempy[2]) options(plot.survfit = list(plotclip=c(tempx, tempy), plotreset=par('usr'))) } @@ -13256,8 +13265,8 @@ if (argclass[ii] == "tdc") { if (argname[[ii]] %in% tevent) stop("attempt to turn event variable", argname[[ii]], "into a tdc") - #if (!is.null(newvar)) - #warning(paste0("replacement of variable '", argname[ii], "'")) + if (!is.null(newvar)) + warning(paste0("replacement of variable '", argname[ii], "'")) # id can be any data type; feed integers to the C routine storage.mode(dstop) <- storage.mode(etime) <- "double" #if time is integer diff -Nru survival-3.1-8/noweb/code.toc survival-3.1-11/noweb/code.toc --- survival-3.1-8/noweb/code.toc 2019-12-01 17:13:38.000000000 +0000 +++ survival-3.1-11/noweb/code.toc 2020-03-05 13:32:34.000000000 +0000 @@ -1,7 +1,7 @@ \contentsline {section}{\numberline {1}Introduction}{2}{section.1} \contentsline {section}{\numberline {2}Cox Models}{2}{section.2} \contentsline {subsection}{\numberline {2.1}Coxph}{2}{subsection.2.1} -\contentsline {subsection}{\numberline {2.2}Exact partial likelihood}{26}{subsection.2.2} +\contentsline {subsection}{\numberline {2.2}Exact partial likelihood}{27}{subsection.2.2} \contentsline {subsection}{\numberline {2.3}Andersen-Gill fits}{36}{subsection.2.3} \contentsline {subsection}{\numberline {2.4}Predicted survival}{55}{subsection.2.4} \contentsline {subsubsection}{\numberline {2.4.1}Multi-state models}{78}{subsubsection.2.4.1} @@ -17,7 +17,7 @@ \contentsline {section}{\numberline {7}Accelerated Failure Time models}{166}{section.7} \contentsline {subsection}{\numberline {7.1}Residuals}{166}{subsection.7.1} \contentsline {section}{\numberline {8}Survival curves}{174}{section.8} -\contentsline {subsection}{\numberline {8.1}Kaplan-Meier}{178}{subsection.8.1} +\contentsline {subsection}{\numberline {8.1}Kaplan-Meier}{179}{subsection.8.1} \contentsline {subsection}{\numberline {8.2}Kaplan-Meier}{183}{subsection.8.2} \contentsline {subsection}{\numberline {8.3}Competing risks}{205}{subsection.8.3} \contentsline {subsubsection}{\numberline {8.3.1}C-code}{215}{subsubsection.8.3.1} diff -Nru survival-3.1-8/noweb/concordance.Rnw survival-3.1-11/noweb/concordance.Rnw --- survival-3.1-8/noweb/concordance.Rnw 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/noweb/concordance.Rnw 2020-01-20 14:14:22.000000000 +0000 @@ -401,7 +401,7 @@ else stop("left hand side of the formula must be a numeric vector or a survival object") } if (timefix) y <- aeqSurv(y) - rval <- list(y= y, x= predict(object)) + rval <- list(y= y, x= predict(object, newdata)) # the type of prediction does not matter, as long as it is a # monotone transform of the linear predictor } @@ -669,7 +669,7 @@ if (is.null(wt)) vmat <- crossprod(dfbeta) else vmat <- t(wt * dfbeta) %*% dfbeta rval <- list(concordance=concordance, count=count, - n=length(flist[[1]]$x), var=vmat, + n=flist[[1]]$n, var=vmat, cvar= sapply(flist, function(x) x$cvar)) if (influence==1) rval$dfbeta <- dfbeta @@ -818,7 +818,7 @@ Referring again to figure \ref{treefig}, \code{btree(13)} yields the vector \code{7 3 8 1 9 4 10 0 11 5 12 2 6} meaning that the smallest element -will be in position 8 of the tree, the next smallest in position 2, etc, +will be in position 8 of the tree, the next smallest in position 4, etc, and using indexing that starts at 0 since the results will be passed to a C routine. The code just above takes care to do all arithmetic as integer. diff -Nru survival-3.1-8/noweb/coxph.Rnw survival-3.1-11/noweb/coxph.Rnw --- survival-3.1-8/noweb/coxph.Rnw 2019-12-01 17:39:58.000000000 +0000 +++ survival-3.1-11/noweb/coxph.Rnw 2020-03-05 13:29:54.000000000 +0000 @@ -323,7 +323,7 @@ Those rows can be removed from the model frame before creating the X matrix. Also identify partially used rows, ones where the necessary covariates are present for some of the possible transitions but not all. -Those obs are dealt with later by the stackup function. +Those obs are dealt with later by the stacker function. <>= # first vector will be true if there is at least 1 transition for which all # covariates are present, second if there is at least 1 for which some are not @@ -370,15 +370,22 @@ Also create the initial values vector. The stacker function will create a separate block of observations for every -unique value in \code{cmat[1,]}. +unique value in \code{cmap[1,]}. Now say that two transitions A:B and A:C share the same baseline hazard. Then either a B or a C outcome will be an ``event'' in that stratum; they would only be distinguished by perhaps having different covariates. +The first thing we do with the result is to rebuild the transitions matrix: +the current version was created before removing missings, and can +seriously overstate the number of transitions available. +Then set up the data. <>= cmap <- parsecovar3(tmap, colnames(X), attr(X, "assign")) xstack <- stacker(cmap, as.integer(istate), X, Y, strata=istrat, states=states) +rkeep <- unique(xstack$rindex) +transitions <- survcheck2(Y[rkeep,], id[rkeep], istate[rkeep])$transitions + X <- xstack$X Y <- xstack$Y istrat <- xstack$strata diff -Nru survival-3.1-8/noweb/coxph.Rnw.orig survival-3.1-11/noweb/coxph.Rnw.orig --- survival-3.1-8/noweb/coxph.Rnw.orig 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/noweb/coxph.Rnw.orig 2020-03-05 13:27:33.000000000 +0000 @@ -64,10 +64,8 @@ tcl <- attr(Terms, 'specials')$cluster if (length(tcl) > 1) stop("a formula cannot have multiple cluster terms") - mycall <- Call + if (length(tcl) > 0) { # there is one - if (!is.null(Call$cluster)) - stop("cannot have both cluster() in the formula and the cluster option") # subscripting of formulas is broken at least through R 3.5, if the # formula contains an offset. Adding offset to the "specials" above # is just a sneaky way to find out if one is present, then call @@ -91,12 +89,15 @@ if (oo < tcl) temp <- c(ooterm, temp) else temp <- c(temp, ooterm) } - mycall$cluster <- attr(Terms, "variables")[[1+tcl]][[2]] + if (is.null(Call$cluster)) + Call$cluster <- attr(Terms, "variables")[[1+tcl]][[2]] + else warning("cluster appears both in a formula and as an argument, formula term ignored") if (is.list(formula)) formula[[1]][[3]] <- reformulate(temp[1-tcl])[[2]] else formula[[3]] <- reformulate(temp[1-tcl])[[2]] - mycall$formula <- formula + Call$formula <- formula + } # create a call to model.frame() that contains the formula (required) @@ -104,9 +105,9 @@ # but don't evaluate it just yet indx <- match(c("formula", "data", "weights", "subset", "na.action", "cluster", "id", "istate"), - names(mycall), nomatch=0) + names(Call), nomatch=0) if (indx[1] ==0) stop("A formula argument is required") - tform <- mycall[c(1,indx)] # only keep the arguments we wanted + tform <- Call[c(1,indx)] # only keep the arguments we wanted tform[[1L]] <- quote(stats::model.frame) # change the function called # if the formula is a list, do the first level of processing on it. @@ -150,12 +151,6 @@ if (!multi && multiform) stop("formula is a list but the response is not multi-state") - # Grab the id variable to check out multi-state data - id <- model.extract(mf, "id") - if (multi) { - <> - } - if (control$timefix) Y <- aeqSurv(Y) <> @@ -164,13 +159,28 @@ # should be extracted after the transform # strats <- attr(Terms, "specials")$strata + hasinteractions <- FALSE + dropterms <- NULL 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) istrat <- as.integer(strata.keep) + + for (i in stemp$vars) { #multiple strata terms are allowed + # The factors attr has one row for each variable in the frame, one + # col for each term in the model. Pick rows for each strata + # var, and find if it participates in any interactions. + if (any(attr(Terms, 'order')[attr(Terms, "factors")[i,] >0] >1)) + hasinteractions <- TRUE + } + 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)) { @@ -194,6 +204,14 @@ } contrast.arg <- NULL #due to shared code with model.matrix.coxph + attr(Terms, "intercept") <- 1 # always have a baseline hazard + + # Grab the id variable to check out multi-state data + id <- model.extract(mf, "id") + if (multi) { + <> + } + <> <> if (multi) { @@ -265,6 +283,12 @@ @ <>= +# remove strata from dformula, before any further processing +if (length(dropterms)){ + Terms2 <- Terms[-dropterms] + dformula <- formula(Terms2) +} else Terms2 <- Terms + # check for consistency of the states, and create a transition # matrix if (length(id)==0) @@ -283,9 +307,9 @@ # build tmap, which has one row per term, one column per transition if (missing(statedata)) covlist2 <- parsecovar2(covlist, NULL, dformula= dformula, - Terms, transitions, states) + Terms2, transitions, states) else covlist2 <- parsecovar2(covlist, statedata, dformula= dformula, - Terms, transitions, states) + Terms2, transitions, states) tmap <- covlist2$tmap if (!is.null(covlist)) { <> @@ -294,7 +318,7 @@ For multi-state models we can't tell what observations should be removed until any extra formulas have been processed. -There may still be rows that are missing \emph{some} of the covariates but +There may be rows that are missing \emph{some} of the covariates but are okay for \emph{some} transitions. Others could be useless. Those rows can be removed from the model frame before creating the X matrix. Also identify partially used rows, ones where the necessary covariates are @@ -799,20 +823,6 @@ We instead assume that nothing is reordered and do a shift. <>= -attr(Terms, "intercept") <- 1 -stemp <- untangle.specials(Terms, 'strata', 1) -hasinteractions <- FALSE -dropterms <- NULL -if (length(stemp$vars) > 0) { #if there is a strata statement - for (i in stemp$vars) { #multiple strata terms are allowed - # The factors attr has one row for each variable in the frame, one - # col for each term in the model. Pick rows for each strata - # var, and find if it participates in any interactions. - if (any(attr(Terms, 'order')[attr(Terms, "factors")[i,] >0] >1)) - hasinteractions <- TRUE - } - if (!hasinteractions) dropterms <- stemp$terms -} if (length(dropterms)) { Terms2 <- Terms[ -dropterms] @@ -865,7 +875,8 @@ <>= if (sum(Y[, ncol(Y)]) == 0) { # No events in the data! - ncoef < ctemp <- rep (NA, ncoef) + 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) @@ -878,7 +889,7 @@ residuals = rep(0.0, data.n), means = colMeans(X), method=method, n = data.n, nevent=0, terms=Terms, assign=assign, - concordance=concordance, + concordance=concordance, wald.test=0.0, y = Y, call=Call) class(rval) <- "coxph" return(rval) @@ -1011,6 +1022,7 @@ else if (length(strats)) fit$strata <- strata.keep } if (y) fit$y <- Y + fit$timefix <- control$timefix # remember this option } @ If any of the weights were not 1, save the results. @@ -1028,6 +1040,7 @@ fit$states <- states fit$cmap <- cmap fit$resid <- rowsum(fit$resid, xstack$rindex) + names(fit$coefficients) <- seq(along=fit$coefficients) if (x) fit$strata <- istrat # save the expanded strata class(fit) <- c("coxphms", class(fit)) } @@ -1070,6 +1083,24 @@ } else dropterms <- NULL + strats <- attr(Terms, "specials")$strata + hasinteractions <- FALSE + 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) + istrat <- as.integer(strata.keep) + + for (i in stemp$vars) { #multiple strata terms are allowed + # The factors attr has one row for each variable in the frame, one + # col for each term in the model. Pick rows for each strata + # var, and find if it participates in any interactions. + if (any(attr(Terms, 'order')[attr(Terms, "factors")[i,] >0] >1)) + hasinteractions <- TRUE + } + if (!hasinteractions) dropterms <- c(dropterms, stemp$terms) + } else istrat <- NULL + <> X } diff -Nru survival-3.1-8/noweb/maximacode survival-3.1-11/noweb/maximacode --- survival-3.1-8/noweb/maximacode 1970-01-01 00:00:00.000000000 +0000 +++ survival-3.1-11/noweb/maximacode 2020-03-05 13:29:54.000000000 +0000 @@ -0,0 +1,102 @@ +/* all of this code was copy-pasted into the maxima program, + and results used to build the R program */ + +/* 2 states, 2 transitions + Remember that R stores matrices columnwise, when copying these results + over. This code uses a,b, c,... in rows. +*/ + +tmat : matrix([-a,a], [b,-b]); +etmat: matrixexp(tmat, t); + +d1: a+b; +e1 : exp(-d1*t); +test: matrix([(a*e1 +b)/d1, (a-a*e1)/d1], [(b-b*e1)/d1, (a+b*e1)/d1]); +expand(test- etmat); /* should be all zeros */ + +deriv1a: diff(etmat[1,1], a); +test: e1*(1-a*t)/d1 - (a*e1 +b)/d1^2; +expand(test - deriv1a); /* should print a 0 */ + +deriv2a: diff(etmat[1,2], a); +deriv3a: diff(etmat[2,1], a); +expand(deriv3a - (b/d1)*(t*e1 + (e1-1)/d1)); + +deriv4a: diff(etmat[2,2], a); + +deriv1b: diff(etmat[1,1], b); +deriv2b: diff(etmat[1,2], b); +deriv3b: diff(etmat[2,1], b); +deriv4b: diff(etmat[2,2], b); + +/* 3 states competing risks */ +tmat : matrix([-(a+b), a, b], [0,0,0], [0,0,0]); +etmat: matrixexp(tmat, t); + +d1: a + b; +e1: exp(-d1*t); +test: matrix([e1, a*(1-e1)/d1, b*(1-e1)/d1], [0,1, 0], [0,0, 1]); +expand(test-etmat); /* should print 0 */ + +dmata: diff(etmat, a); +dmatb: diff(etmat, b); + +test: (1+ a*t*e1 - e1)/d1 + a(e1-1)/d1^2; +expand(test - dmata[1,2]); + +test: (b/d1)*(t*e1 + (e1-1)/d1); +expand(test - dmata[1,3]); + +test: (a/d1)*(t*e1 + (e1-1)/d1); +expand(test - dmatb[1,2]); + +/* 4 states, competing risks */ +tmat : matrix([-(a+b+d), a, b, d], [0,0,0,0], [0,0,0,0], [0,0,0,0]); +etmat: matrixexp(tmat, t); +dmata: diff(etmat, a); +dmatb: diff(etmat, b); + +d1: a+b+d; +e1: exp(-d1*t); +test: (1+ a*t*e1 - e1)/d1 + a(e1-1)/d1^2; +expand(test - dmata[1,2]); + +test: (b/d1)*(t*e1 + (e1-1)/d1); +expand(test - dmata[1,3]); + + +/* 3 states, a to b and all to death */ +tmat : matrix([-(a+b), a, b], [0, -c, c], [0,0,0]); +etmat: matrixexp(tmat * t); + +d1 : c- a+b; +e1 : exp(-(a+b)*t); +e2 : exp(-c*t); +expand(a*(e1 - e2)/d1 - etmat[1,2]); /* simpler form for the term */ +expand((b-c)*e1/d1 + a*e2/d1 +1 - etmat[1,3]); + +expand(diff(etmat[1,1], a)); +diff(etmat[1,2], a); +diff(etmat[2,2], a); + +test3: (e1-e2)*(1+ a/d2)/d2 - a*t*e1/d2; +expand(test3- diff(etmat[1,2], a)); +test2: t + a*t/d2 - +diff(etmat[2,2], a) - (test1 + test2); + +expand(diff(etmat[1,3], a)); +expand(diff(etmat[2,2], a)); + + +/* special case of the above, where c= a+b */ +tmat : matrix([-(a+b), a, b], [0,-(a+b), a+b], [0,0,0]); +etmat: matrixexp(tmat*t); + +expand(diff(etmat[1,1], a)); +expand(diff(etmat[1,2], a)); +expand(diff(etmat[1,3], a)); +expand(diff(etmat[2,2], a)) + +/* illness-death */ +tmat : matrix([-(a+b), a, b], [c, -(c+d), d], [0,0,0]); +etmat: matrixexp(tmat*t); diff -Nru survival-3.1-8/noweb/msurv.nw survival-3.1-11/noweb/msurv.nw --- survival-3.1-8/noweb/msurv.nw 2019-12-01 16:48:18.000000000 +0000 +++ survival-3.1-11/noweb/msurv.nw 2020-02-03 19:28:02.000000000 +0000 @@ -1,6 +1,6 @@ \subsubsection{C-code} (This is set up as a separate file in the source code directory since -it is easier to make the code stay in C-mode if the file has a .nw +it is easier to make emacs stay in C-mode if the file has a .nw extension.) <>= diff -Nru survival-3.1-8/noweb/plot.Rnw survival-3.1-11/noweb/plot.Rnw --- survival-3.1-8/noweb/plot.Rnw 2019-10-31 13:16:45.000000000 +0000 +++ survival-3.1-11/noweb/plot.Rnw 2020-01-30 17:51:40.000000000 +0000 @@ -197,7 +197,7 @@ else stop("invalid cumhaz argument") } else if (inherits(x, "survfitms")) { - i <- (x$states != noplot) + i <- !(x$states %in% noplot) if (all(i) || !any(i)) { # the !any is a failsafe, in case none are kept then ignore noplot ssurv <- smat(x$pstate) @@ -383,7 +383,7 @@ if (!missing(xlim) && !is.null(xlim)) { tempx <- xlim - if (xaxs == 'S') tempx[2] <- tempx[2] + diff(tempx)*1.04 + if (xaxs == 'S') tempx[2] <- tempx[1] + diff(tempx)*1.04 } else { temp <- stime[is.finite(stime)] @@ -454,7 +454,7 @@ if (resetclip) { # yes, do it - if (xaxs=='S') tempx <- c(tempx[1], temp[1]) + #if (xaxs=='S') tempx <- c(tempx[1], temp[1]) clip(tempx[1], tempx[2], tempy[1], tempy[2]) options(plot.survfit = list(plotclip=c(tempx, tempy), plotreset=par('usr'))) } diff -Nru survival-3.1-8/noweb/residuals.survreg.Rnw survival-3.1-11/noweb/residuals.survreg.Rnw --- survival-3.1-8/noweb/residuals.survreg.Rnw 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/noweb/residuals.survreg.Rnw 2020-02-21 19:41:52.000000000 +0000 @@ -340,7 +340,9 @@ else score <- cbind(score, deriv[,4]) } rr <- score %*% vv - if (type=='dfbetas') rr <- rr %*% diag(1/sqrt(diag(vv))) + # cause column names to be retained + # old: if (type=='dfbetas') rr[] <- rr %*% diag(1/sqrt(diag(vv))) + if (type=='dfbetas') rr <- rr * rep(1/sqrt(diag(vv)), each=nrow(rr)) if (type=='ldcase') rr<- rowSums(rr*score) } diff -Nru survival-3.1-8/noweb/tmerge.Rnw survival-3.1-11/noweb/tmerge.Rnw --- survival-3.1-8/noweb/tmerge.Rnw 2019-11-20 19:27:36.000000000 +0000 +++ survival-3.1-11/noweb/tmerge.Rnw 2020-01-28 19:06:11.000000000 +0000 @@ -479,8 +479,8 @@ if (argclass[ii] == "tdc") { if (argname[[ii]] %in% tevent) stop("attempt to turn event variable", argname[[ii]], "into a tdc") - #if (!is.null(newvar)) - #warning(paste0("replacement of variable '", argname[ii], "'")) + if (!is.null(newvar)) + warning(paste0("replacement of variable '", argname[ii], "'")) # id can be any data type; feed integers to the C routine storage.mode(dstop) <- storage.mode(etime) <- "double" #if time is integer diff -Nru survival-3.1-8/noweb/tmerge.Rnw.orig survival-3.1-11/noweb/tmerge.Rnw.orig --- survival-3.1-8/noweb/tmerge.Rnw.orig 1970-01-01 00:00:00.000000000 +0000 +++ survival-3.1-11/noweb/tmerge.Rnw.orig 2020-01-28 15:18:01.000000000 +0000 @@ -0,0 +1,626 @@ +\section{tmerge} +The tmerge function was designed around a set of specific problems. +The idea is to build up a time dependent data set one endpoint at at time. +The primary arguments are +\begin{itemize} + \item data1: the base data set that will be added onto + \item data2: the source for new information + \item id: the subject identifier in the new data + \item \ldots: additional arguments that add variables to the data set + \item tstart, tstop: used to set the time range for each subject + \item options +\end{itemize} +The created data set has three new variables (at least), which are +\code{id}, \code{tstart} and \code{tstop}. + +The key part of the call are the ``\ldots'' arguments +which each can be one of four types: +tdc() and cumtdc() add a time dependent variable, event() and cumevent() +add a new endpoint. +In the survival routines time intervals are open on the left and +closed on the right, i.e., (tstart, tstop]. +Time dependent covariates apply from the start of an interval and events +occur at the end of an interval. +If a data set already had intervals of (0,10] and (10, 14] a new time +dependent covariate or event at time 8 would lead to three intervals of +(0,8], (8,10], and (10,14]; +the new time-dependent covariate value would be added to the second interval, +a new event would be added to the first one. + +A typical call would be +<>= + newdata <- tmerge(newdata, old, id=clinic, diabetes=tdc(diab.time)) +@ +which would add a new time dependent covariate \code{diabetes} to the +data set. + +<>= +tmerge <- function(data1, data2, id, ..., tstart, tstop, options) { + Call <- match.call() + # The function wants to recognize special keywords in the + # arguments, so define a set of functions which will be used to + # mark objects + new <- new.env(parent=parent.frame()) + assign("tdc", function(time, value=NULL) { + x <- list(time=time, value=value); + class(x) <- "tdc"; x}, + envir=new) + assign("cumtdc", function(time, value=NULL) { + x <- list(time=time, value=value); + class(x) <-"cumtdc"; x}, + envir=new) + assign("event", function(time, value=NULL, censor=NULL) { + x <- list(time=time, value=value, censor=censor); + class(x) <-"event"; x}, + envir=new) + assign("cumevent", function(time, value=NULL, censor=NULL) { + x <- list(time=time, value=value, censor=censor); + class(x) <-"cumevent"; x}, + envir=new) + + if (missing(data1) || missing(data2) || missing(id)) + stop("the data1, data2, and id arguments are required") + if (!inherits(data1, "data.frame")) stop("data1 must be a data frame") + + <> + <> + <> +} +@ + +The program can't use formulas because the \ldots arguments need to be +named. This results in a bit of evaluation magic to correctly assess +arguments. +The routine below could have been set out as a separate top-level routine, +the argument is where we want to document it: within the tmerge page or not. +I decided on the former. +<>= +tmerge.control <- function(idname="id", tstartname="tstart", tstopname="tstop", + delay =0, na.rm=TRUE, tdcstart=NA, ...) { + extras <- list(...) + if (length(extras) > 0) + stop("unrecognized option(s):", paste(names(extras), collapse=', ')) + if (length(idname) != 1 || make.names(idname) != idname) + stop("idname option must be a valid variable name") + if (!is.null(tstartname) && + (length(tstartname) !=1 || make.names(tstartname) != tstartname)) + stop("tstart option must be NULL or a valid variable name") + if (length(tstopname) != 1 || make.names(tstopname) != tstopname) + stop("tstop option must be a valid variable name") + if (length(delay) !=1 || !is.numeric(delay) || delay < 0) + stop("delay option must be a number >= 0") + if (length(na.rm) !=1 || ! is.logical(na.rm)) + stop("na.rm option must be TRUE or FALSE") + if (length(tdcstart) !=1) stop("tdcstart must be a single value") + list(idname=idname, tstartname=tstartname, tstopname=tstopname, + delay=delay, na.rm=na.rm, tdcstart=tdcstart) +} + +tname <- attr(data1, "tname") +firstcall <- is.null(tname) #first call to the function +if (!firstcall && any(is.null(match(unlist(tname), names(data1))))) + stop("data1 does not match its own tname attribute") + +if (!missing(options)) { + if (!is.list(options)) stop("options must be a list") + if (!is.null(tname)) { + # If an option name matches one already in tname, don't confuse + # the tmerge.control routine with duplicate arguments + temp <- match(names(options), names(tname), nomatch=0) + topt <- do.call(tmerge.control, c(options, tname[temp==0])) + if (any(temp >0)) { + # A variable name is changing midstream, update the + # variable names in data1 + varname <- tname[c("idname", "tstartname", "tstopname")] + temp2 <- match(varname, names(data1)) + names(data1)[temp2] <- varname + } + } + else topt <- do.call(tmerge.control, options) +} +else if (length(tname)) topt <- do.call(tmerge.control, tname) +else topt <- tmerge.control() + +# id, tstart, tstop are found in data2 +if (missing(id)) stop("the id argument is required") +if (missing(data1) || missing(data2)) + stop("two data sets are required") +id <- eval(Call[["id"]], data2, enclos=emptyenv()) #don't find it elsewhere +if (is.null(id)) stop("id variable not found in data2") +if (any(is.na(id))) stop("id variable cannot have missing values") + +if (firstcall) { + if (!missing(tstop)) { + tstop <- eval(Call[["tstop"]], data2) + if (length(tstop) != length(id)) + stop("tstop and id must be the same length") + # The neardate routine will check for legal tstop data type + } + if (!missing(tstart)) { + tstart <- eval(Call[["tstart"]], data2) + if (length(tstart)==1) tstart <- rep(tstart, length(id)) + if (length(tstart) != length(id)) + stop("tstart and id must be the same length") + if (any(tstart >= tstop)) + stop("tstart must be < tstop") + } +} +else { + if (!missing(tstart) || !missing(tstop)) + stop("tstart and tstop arguments only apply to the first call") +} +@ + +Get the \ldots arguments. They are evaluated in a special frame, +set up earlier, so that the definitions of the functions tdc, +cumtdc, event, and cumevent are local to tmerge. +Check that they are all legal: each argument is named, and is one of the four +allowed types. +<>= +# grab the... arguments +notdot <- c("data1", "data2", "id", "tstart", "tstop", "options") +dotarg <- Call[is.na(match(names(Call), notdot))] +dotarg[[1]] <- as.name("list") # The as-yet dotarg arguments +if (missing(data2)) args <- eval(dotarg, envir=new) +else args <- eval(dotarg, data2, enclos=new) + +argclass <- sapply(args, function(x) (class(x))[1]) +argname <- names(args) +if (any(argname== "")) stop("all additional argments must have a name") + +check <- match(argclass, c("tdc", "cumtdc", "event", "cumevent")) +if (any(is.na(check))) + stop(paste("argument(s)", argname[is.na(check)], + "not a recognized type")) +@ + +The tcount matrix keeps track of what we have done, and is added to +the final object at the end. +This is useful to the user for debugging what may have gone right or +wrong in their usage. + +<>= +# The tcount matrix is useful for debugging +tcount <- matrix(0L, length(argname), 8) +dimnames(tcount) <- list(argname, c("early","late", "gap", "within", + "boundary", "leading", "trailing", + "tied")) +tevent <- attr(data1, "tevent") # event type variables +tcens <- attr(data1, "tcensor")# censor code for variables +if (is.null(tcens)) tcens <- vector('list', 0) +@ + +The very first call to the routine is special, since this is when the +range of legal times is set. We also apply an initial sort to the data +if necessary so that times are in order. +There are 2 cases: +\begin{enumerate} + \item Adding a time range: tstop comes from data2, optional tstart, and the + id can be simply matched, by which we mean no duplicates in data1. + \item The more common case: there is no tstop, one observation per subject, + and the first optional argument is an + event or cumevent. We then use its time as the range. +\end{enumerate} +One thing we could add, but didn't, was to warn if any of the three new +variables will stomp on ones already in data1. + +<>= +newdata <- data1 #make a copy +if (firstcall) { + # We don't look for topt$id. What if the user had id=clinic, but their + # starting data set also had a variable named "id". We want clinic for + # this first call. + idname <- Call[["id"]] + if (!is.name(idname)) + stop("on the first call 'id' must be a single variable name") + + # The line below finds tstop and tstart variables in data1 + indx <- match(c(topt$idname, topt$tstartname, topt$tstopname), names(data1), + nomatch=0) + if (any(indx[1:2]>0) && FALSE) { # warning currently turned off. Be chatty? + overwrite <- c(topt$tstartname, topt$tstopname)[indx[2:3]] + warning("overwriting data1 variables", paste(overwrite, collapse=' ')) + } + + temp <- as.character(idname) + if (!is.na(match(temp, names(data1)))) { + data1[[topt$idname]] <- data1[[temp]] + baseid <- data1[[temp]] + } + else stop("id variable not found in data1") + + if (any(duplicated(baseid))) + stop("for the first call (that establishes the time range) data1 must have no duplicate identifiers") + + if (length(baseid)== length(id) && all(baseid == id)) newdata <- data1 + else { # Note: 'id' is the idlist for data 2 + indx2 <- match(id, baseid) + if (any(is.na(indx2))) + stop("'id' has values not in data1") + newdata <- data1[indx2,] + } + if (missing(tstop)) { # case 2 + if (length(argclass)==0 || argclass[1] != "event") + stop("neither a tstop argument nor an initial event argument was found") + tstop <- args[[1]][[1]] + } + + # at this point newdata and data2 are in the same order, same # rows + if (any(is.na(tstop))) + stop("missing time value, when that variable defines the span") + if (missing(tstart)) { + indx <- which(tstop <=0) + if (length(indx) >0) stop("found an ending time of ", tstop[indx[1]], + ", the default starting time of 0 is invalid") + tstart <- rep(0, length(tstop)) + } + if (any(tstart >= tstop)) + stop("tstart must be < tstop") + newdata[[topt$tstartname]] <- tstart + newdata[[topt$tstopname]] <- tstop + n <- nrow(newdata) + if (any(duplicated(id))) { + # sort by time within id + indx1 <- match(id, unique(id)) + newdata <- newdata[order(indx1, tstop),] + } + temp <- newdata[[topt$idname]] + if (any(tstart >= tstop)) stop("tstart must be < tstop") + if (any(newdata$tstart[-n] > newdata$tstop[-1] & + temp[-n] == temp[-1])) + stop("first call has created overlapping or duplicated time intervals") +} +else { #not a first call + if (any(is.na(match(id, data1[[topt$idname]])))) + stop("id values were found in data2 which are not in data1") +} +@ + +Now for the real work. For each additional argument we first match the +id/time pairs of the new data to the current data set, and categorize +each into a type. If the time value in data2 is NA, then that +addition is skipped. Ditto if the value is NA and options narm=TRUE. +This is a convenience for the user, who will often +be merging in a variable like ``day of first diabetes diagnosis'' which +is missing for those who never had that outcome occur. +<>= +saveid <- id +for (ii in seq(along.with=args)) { + argi <- args[[ii]] + baseid <- newdata[[topt$idname]] + dstart <- newdata[[topt$tstartname]] + dstop <- newdata[[topt$tstopname]] + argcen <- argi$censor + + # if an event time is missing then skip that obs + etime <- argi$time + if (length(etime) != length(saveid)) + stop("argument ", argname[ii], " is not the same length as id") + if (!is.null(argi$value)) { + if (length(argi$value) != length(saveid)) + stop("argument ", argname[ii], " is not the same length as id") + if (topt$na.rm) keep <- !(is.na(etime) | is.na(argi$value)) + else keep <- !is.na(etime) + if (!all(keep)) { + etime <- etime[keep] + argi$value <- argi$value[keep] + } + } + else { + keep <- !is.na(etime) + etime <- etime[keep] + } + id <- saveid[keep] + + # Later steps become easier if we sort the new data by id and time + # The match() is critical when baseid is not in sorted order. The + # etime part of the sort will change from one ii value to the next. + indx <- order(match(id, baseid), etime) + id <- id[indx] + etime <- etime[indx] + if (!is.null(argi$value)) + yinc <- argi$value[indx] + else yinc <- NULL + + # indx1 points to the closest start time in the baseline data (data1) + # that is <= etime. indx2 to the closest end time that is >=etime. + # If etime falls into a (tstart, tstop) interval, indx1 and indx2 + # will match + # If the "delay" argument is set and this event is of type tdc, then + # move any etime that is after the entry time for a subject. + if (topt$delay >0 && argclass[ii] %in% c("tdc", "cumtdc")) { + mintime <- tapply(dstart, baseid, min) + index <- match(id, names(mintime)) + etime <- ifelse(etime <= mintime[index], etime, etime+ topt$delay) + } + + indx1 <- neardate(id, baseid, etime, dstart, best="prior") + indx2 <- neardate(id, baseid, etime, dstop, best="after") + + # The event times fall into one of 5 categories + # 1. Before the first interval + # 2. After the last interval + # 3. Outside any interval but with time span, i.e, it falls into + # a gap in follow-up + # 4. Strictly inside an interval (does't touch either end) + # 5. Inside an interval, but touching. + itype <- ifelse(is.na(indx1), 1, + ifelse(is.na(indx2), 2, + ifelse(indx2 > indx1, 3, + ifelse(etime== dstart[indx1] | + etime== dstop[indx2], 5, 4)))) + + # Subdivide the events that touch on a boundary + # 1: intervals of (a,b] (b,d], new count at b "tied edge" + # 2: intervals of (a,b] (c,d] with c>b, new count at c, "front edge" + # 3: intervals of (a,b] (c,d] with c>b, new count at b, "back edge" + # + subtype <- ifelse(itype!=5, 0, + ifelse(indx1 == indx2+1, 1, + ifelse(etime==dstart[indx1], 2, 3))) + tcount[ii,1:7] <- table(factor(itype+subtype, levels=c(1:4, 6:8))) + + # count ties. id and etime are not necessarily sorted + tcount[ii,8] <- sum(tapply(etime, id, function(x) sum(duplicated(x)))) + + <> +} +@ + +A \code{tdc} or \code{cumtdc} operator defines a new time-dependent +variable which applies to all future times. +Say that we had the following scenario for one subject +\begin{center} + \begin{tabular}{rr|rr} + \multicolumn{2}{c}{current} & \multicolumn{2}{c}{addition} \\ + tstart & tstop & time & x \\ + 2 & 5 & 1 & 20.2 \\ + 6 & 7 & 7 & 11 \\ + 7 & 15 & 8 & 17.3 \\ + 15 & 30 \\ + \end{tabular} + \end{center} +The resulting data set will have intervals of (2,5), (6,7), (7,8) and (8,15) +with covariate values of 20.2, 20.2, 11, and 17.3. +Only a covariate change that occurs within an interval causes a new data +row. Covariate changes that happen after the last interval are ignored, +i.e. at change at time $\ge 30$ in the above example. + +If instead this had been events at times 1, 7, and 8, the first event would +be ignored since it happens outside of any interval, so would an event +at exactly time 2. The event at time +7 would be recorded in the (6,7) interval and the one at time 8 in the +(7,8) interval: events happen at the ends of intervals. +In both cases new rows are only generated for new time values that fall +strictly within one of the old intervals. + +When a subject has two increments on the same day the later one wins. +This is correct behavior for cumtdc, a bit odd for cumevent, and the +user's problem for tdc and event. +We report back the number of ties so that the user can deal with it. + +Where are we now with the variables? +\begin{center} + \begin{tabular}{cccc} + itype& class & indx1 & indx2 \\ \hline + 1 & before & NA & next interval \\ + 2 & after & prior interval & NA \\ + 3 & in a gap & prior interval & next interval \\ + 4 & within interval & containing interval & containing interval \\ + 5-1 & on a join & next interval & prior interval \\ + 5-2 & front edge & containing & containing \\ + 5-3 & back edge & containing & containing \\ + \end{tabular} +\end{center} +If there are any itype 4, start by expanding the data set to add +new cut points, which will turn all the 4's into 5-1 types. +When expanding, all the event type variables turn into ``censor'' at the +newly added times and other variables stay the same. +A subject could have more than one new cutpoint added within an interval +so we have to count each. +In newdata all the rows for a given subject are contiguous and in time +order, though the data set may not be in subject order. +<>= +indx4 <- which(itype==4) +n4 <- length(indx4) +if (n4 > 0) { + # we need to eliminate duplicate times within the same id, but + # do so without changing the class of etime: it might + # be a Date, an integer, a double, ... + # Using unique on a data.frame does the trick + icount <- data.frame(irow= indx1[indx4], etime=etime[indx4]) + icount <- unique(icount) + # the icount data frame will be sorted by second column within first + # so rle is faster than table + n.add <- rle(icount$irow)$length # number of rows to add for each id + + # expand the data + irep <- rep.int(1L, nrow(newdata)) + erow <- unique(indx1[indx4]) # which rows in newdata to be expanded + irep[erow] <- 1+ n.add # number of rows in new data + jrep <- rep(1:nrow(newdata), irep) #stutter the duplicated rows + newdata <- newdata[jrep,] #expand it out + dstart <- dstart[jrep] + dstop <- dstop[jrep] + + #fix up times + nfix <- length(erow) + temp <- vector("list", nfix) + iend <- (cumsum(irep))[irep >1] #end row of each duplication set + for (j in 1:nfix) temp[[j]] <- -(seq(n.add[j] -1, 0)) + iend[j] + newrows <- unlist(temp) + + dstart[newrows] <- dstop[newrows-1] <- icount$etime + newdata[[topt$tstartname]] <- dstart + newdata[[topt$tstopname]] <- dstop + for (ename in tevent) newdata[newrows-1, ename] <- tcens[[ename]] + + # refresh indices + baseid <- newdata[[topt$idname]] + indx1 <- neardate(id, baseid, etime, dstart, best="prior") + indx2 <- neardate(id, baseid, etime, dstop, best="after") + subtype[itype==4] <- 1 #all the "insides" are now on a tied edge + itype[itype==4] <- 5 +} +@ + +Now we can add the new variable. +The most common is a tdc, so start with it. +The C routine returns a set of indices: 0,1,1,2,3,0,4,... would mean that +row 1 of the new data happens before the tdc variable, 2 and 3 take values from +the first element of yinc, etc. +By returning an index, the yinc variable can be of any data type. Using +is.na() on the left side below causes the \emph{right} kind of NA to be inserted +(this trick was stolen from the merge routine). + +<>= +# add a tdc variable +newvar <- newdata[[argname[ii]]] # prior value (for sequential tmerge calls) +if (argclass[ii] == "tdc") { + if (argname[[ii]] %in% tevent) + stop("attempt to turn event variable", argname[[ii]], "into a tdc") + #if (!is.null(newvar)) + #warning(paste0("replacement of variable '", argname[ii], "'")) + + # id can be any data type; feed integers to the C routine + storage.mode(dstop) <- storage.mode(etime) <- "double" #if time is integer + uid <- unique(baseid) + index <- .Call(Ctmerge2, match(baseid, uid), dstop, + match(id, uid), etime) + + if (is.null(yinc)) { + # add a 0/1 variable + if (is.null(newvar))newvar <- ifelse(index==0, 0L, 1L) + else newvar <- ifelse(index==0, newvar, 1L) + } + else { + if (is.null(newvar)) { + newvar <- yinc[pmax(1L, index)] + if (any(index==0)) { + if (is.na(topt$tdcstart)) is.na(newvar) <- (index==0L) + else { + if (is.numeric(newvar)) + newvar[index=0L] <- as.numeric(topt$tdcstart) + else { + if (is.factor(newvar)) { + # special case: if tdcstart isn't in the set, + # add it to the levels + if (is.na(match(topt$tdcstart, levels(newvar)))) + levels(newvar) <- c(levels(newvar), topt$tdcstart) + } + newvar[index== 0L] <- topt$tdcstart + } + } + } + } + else newvar <- ifelse(index==0, newvar, yinc[pmax(1L, index)]) + } +} +@ + +Events and cumevents are easy because +each affects only one interval. +<>= +# add events +if (argclass[ii] %in% c("cumtdc", "cumevent")) { + if (is.null(yinc)) yinc <- rep(1L, length(id)) + else if (is.logical(yinc)) yinc <- as.numeric(yinc) # allow cumulative T/F + if (!is.numeric(yinc)) stop("invalid increment for cumtdc or cumevent") + ykeep <- (yinc !=0) # ignore the addition of a censoring event + yinc <- unlist(tapply(yinc, match(id, baseid), cumsum)) +} + +if (argclass[ii] %in% c("event", "cumevent")) { + if (!is.null(newvar)) { + if (!argname[ii] %in% tevent) { + #warning(paste0("non-event variable '", argname[ii], "' replaced by an event variable")) + newvar <- NULL + } + else if (!is.null(yinc)) { + if (class(newvar) != class(yinc)) + stop("attempt to update an event variable with a different type") + if (is.factor(newvar) && !all(levels(yinc) %in% levels(newvar))) + stop("attemp to update an event variable and levels do not match") + } + } + + if (is.null(yinc)) yinc <- rep(1L, length(id)) + if (is.null(newvar)) { + if (is.numeric(yinc)) newvar <- rep(0L, nrow(newdata)) + else if (is.factor(yinc)) + newvar <- factor(rep(levels(yinc)[1], nrow(newdata)), + levels(yinc)) + else if (is.character(yinc)) newvar <- rep('', nrow(newdata)) + else if (is.logical(yinc)) newvar <- rep(FALSE, nrow(newdata)) + else stop("invalid value for a status variable") + } + + keep <- (subtype==1 | subtype==3) # all other events are thrown away + if (argclass[ii] == "cumevent") keep <- (keep & ykeep) + newvar[indx2[keep]] <- yinc[keep] + + # add this into our list of 'this is an event type variable' + if (!(argname[ii] %in% tevent)) { + tevent <- c(tevent, argname[[ii]]) + if (is.factor(yinc)) tcens <- c(tcens, list(levels(yinc)[1])) + else if (is.logical(yinc)) tcens <- c(tcens, list(FALSE)) + else if (is.character(yinc)) tcens <- c(tcens, list("")) + else tcens <- c(tcens, list(0)) + names(tcens) <- tevent + } +} + +else if (argclass[ii] == "cumtdc") { # process a cumtdc variable + # I don't have a good way to catch the reverse of this user error + if (argname[[ii]] %in% tevent) + stop("attempt to turn event variable", argname[[ii]], "into a cumtdc") + + keep <- itype != 2 # changes after the last interval are ignored + indx <- ifelse(subtype==1, indx1, + ifelse(subtype==3, indx2+1L, indx2)) + + # we want to pass the right kind of NA to the C code + if (is.na(topt$tdcstart)) topt$tdcstart <- as.numeric(topt$tdcstart) + if (is.null(newvar)) { # not overwriting a prior value + if (is.null(argi$value)) newvar <- rep(0.0, nrow(newdata)) + else newvar <- rep(topt$tdcstart, nrow(newdata)) + } + + # the increment must be numeric + if (!is.numeric(newvar)) + stop("data and options$tdcstart do not agree on data type") + # id can be any data type; feed integers to the C routine + storage.mode(yinc) <- storage.mode(dstop) <- "double" + storage.mode(newvar) <- storage.mode(etime) <- "double" + newvar <- .Call(Ctmerge, match(baseid, baseid), dstop, newvar, + match(id, baseid)[keep], etime[keep], + yinc[keep], indx[keep]) +} + +newdata[[argname[ii]]] <- newvar +@ + +Finish up by adding the attributes and the class +<>= +attr(newdata, "tname") <- topt[c("idname", "tstartname", "tstopname")] +attr(newdata, "tcount") <- rbind(attr(data1, "tcount"), tcount) +if (length(tevent)) { + attr(newdata, "tevent") <- tevent + attr(newdata, "tcensor" ) <- tcens + } +row.names(newdata) <- NULL #These are a mess; kill them off. +# Not that it works: R just assigns new row names. +class(newdata) <- c("data.frame") +newdata +@ + +The print routine is for checking: it simply prints out the attributes. +<>= +print.tmerge <- function(x, ...) { + print(attr(x, "tcount")) +} +"[.tmerge" <- function(x, ..., drop=TRUE){ + class(x) <- "data.frame" + NextMethod(,x) + } +@ diff -Nru survival-3.1-8/R/concordance.R survival-3.1-11/R/concordance.R --- survival-3.1-8/R/concordance.R 2019-12-01 17:13:40.000000000 +0000 +++ survival-3.1-11/R/concordance.R 2020-03-05 13:32:36.000000000 +0000 @@ -350,7 +350,7 @@ else stop("left hand side of the formula must be a numeric vector or a survival object") } if (timefix) y <- aeqSurv(y) - rval <- list(y= y, x= predict(object)) + rval <- list(y= y, x= predict(object, newdata)) # the type of prediction does not matter, as long as it is a # monotone transform of the linear predictor } @@ -572,7 +572,7 @@ if (is.null(wt)) vmat <- crossprod(dfbeta) else vmat <- t(wt * dfbeta) %*% dfbeta rval <- list(concordance=concordance, count=count, - n=length(flist[[1]]$x), var=vmat, + n=flist[[1]]$n, var=vmat, cvar= sapply(flist, function(x) x$cvar)) if (influence==1) rval$dfbeta <- dfbeta diff -Nru survival-3.1-8/R/coxph.R survival-3.1-11/R/coxph.R --- survival-3.1-8/R/coxph.R 2019-12-01 17:40:18.000000000 +0000 +++ survival-3.1-11/R/coxph.R 2020-03-05 13:32:36.000000000 +0000 @@ -413,6 +413,9 @@ cmap <- parsecovar3(tmap, colnames(X), attr(X, "assign")) xstack <- stacker(cmap, as.integer(istate), X, Y, strata=istrat, states=states) + rkeep <- unique(xstack$rindex) + transitions <- survcheck2(Y[rkeep,], id[rkeep], istate[rkeep])$transitions + X <- xstack$X Y <- xstack$Y istrat <- xstack$strata diff -Nru survival-3.1-8/R/plot.cox.zph.R survival-3.1-11/R/plot.cox.zph.R --- survival-3.1-8/R/plot.cox.zph.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/R/plot.cox.zph.R 2019-12-09 15:45:53.000000000 +0000 @@ -3,24 +3,14 @@ ...) { xx <- x$x yy <- x$y -# d <- nrow(yy) - df <- max(df) #error proofing + df <- max(df) # in case df is a vector nvar <- ncol(yy) pred.x <- seq(from=min(xx), to=max(xx), length=nsmo) temp <- c(pred.x, xx) lmat <- ns(temp, df=df, intercept=TRUE) pmat <- lmat[1:nsmo,] # for prediction xmat <- lmat[-(1:nsmo),] - qmat <- qr(xmat) - if (qmat$rank < df) - stop("Spline fit is singular, try a smaller degrees of freedom") - - if (se) { - bk <- backsolve(qmat$qr[1:df, 1:df], diag(df)) - xtx <- bk %*% t(bk) -# seval <- d*((pmat%*% xtx) *pmat) %*% rep(1, df) - seval <- ((pmat%*% xtx) *pmat) %*% rep(1, df) - } + if (missing(ylab)) ylab <- paste("Beta(t) for", dimnames(yy)[[2]]) if (missing(var)) var <- 1:nvar @@ -28,7 +18,7 @@ if (is.character(var)) var <- match(var, dimnames(yy)[[2]]) if (any(is.na(var)) || max(var)>nvar || min(var) <1) stop("Invalid variable requested") - } + } # # Figure out a 'good' set of x-axis labels. Find 8 equally spaced @@ -54,12 +44,33 @@ col <- rep(col, length=2) lwd <- rep(lwd, length=2) lty <- rep(lty, length=2) + + # Now, finally do the work for (i in var) { + # Since release 3.1-6, yy can have missing values. If a covariate is + # constant within a stratum then it's Shoenfeld residual is identially + # zero for all observations in that stratum. These "structural zeros" + # are marked with an NA. They contain no information and should not + # by plotted. Thus we need to do the spline fit one stratum at a time. y <- yy[,i] + keep <- !is.na(y) + if (!all(keep)) y <- y[keep] + + qmat <- qr(xmat[keep,]) + if (qmat$rank < df) { + warning("spline fit is singular, variable ", i, " skipped") + next + } + yhat <- pmat %*% qr.coef(qmat, y) if (resid) yr <-range(yhat, y) else yr <-range(yhat) + if (se) { + bk <- backsolve(qmat$qr[1:df, 1:df], diag(df)) + xtx <- bk %*% t(bk) + seval <- ((pmat%*% xtx) *pmat) %*% rep(1, df) + temp <- 2* sqrt(x$var[i,i]*seval) yup <- yhat + temp ylow<- yhat - temp @@ -69,16 +80,16 @@ if (x$transform=='identity') plot(range(xx), yr, type='n', xlab=xlab, ylab=ylab[i], ...) else if (x$transform=='log') - plot(range(xx), yr, type='n', xlab=xlab, ylab=ylab[i], log='x', - ...) + plot(range(xx[keep]), yr, type='n', xlab=xlab, ylab=ylab[i], + log='x', ...) else { - plot(range(xx), yr, type='n', xlab=xlab, ylab=ylab[i], + plot(range(xx[keep]), yr, type='n', xlab=xlab, ylab=ylab[i], axes=FALSE,...) axis(1, xaxisval, xaxislab) axis(2) box() } - if (resid) points(xx, y) + if (resid) points(xx[keep], y) lines(pred.x, yhat, lty=lty[1], col=col[1], lwd=lwd[1]) if (se) { diff -Nru survival-3.1-8/R/plot.survfit.R survival-3.1-11/R/plot.survfit.R --- survival-3.1-8/R/plot.survfit.R 2019-12-01 17:13:43.000000000 +0000 +++ survival-3.1-11/R/plot.survfit.R 2020-03-05 13:32:39.000000000 +0000 @@ -97,7 +97,7 @@ else stop("invalid cumhaz argument") } else if (inherits(x, "survfitms")) { - i <- (x$states != noplot) + i <- !(x$states %in% noplot) if (all(i) || !any(i)) { # the !any is a failsafe, in case none are kept then ignore noplot ssurv <- smat(x$pstate) @@ -258,7 +258,7 @@ if (!missing(xlim) && !is.null(xlim)) { tempx <- xlim - if (xaxs == 'S') tempx[2] <- tempx[2] + diff(tempx)*1.04 + if (xaxs == 'S') tempx[2] <- tempx[1] + diff(tempx)*1.04 } else { temp <- stime[is.finite(stime)] @@ -320,7 +320,7 @@ if (resetclip) { # yes, do it - if (xaxs=='S') tempx <- c(tempx[1], temp[1]) + #if (xaxs=='S') tempx <- c(tempx[1], temp[1]) clip(tempx[1], tempx[2], tempy[1], tempy[2]) options(plot.survfit = list(plotclip=c(tempx, tempy), plotreset=par('usr'))) } @@ -527,7 +527,7 @@ else stop("invalid cumhaz argument") } else if (inherits(x, "survfitms")) { - i <- (x$states != noplot) + i <- !(x$states %in% noplot) if (all(i) || !any(i)) { # the !any is a failsafe, in case none are kept then ignore noplot ssurv <- smat(x$pstate) @@ -871,7 +871,7 @@ else stop("invalid cumhaz argument") } else if (inherits(x, "survfitms")) { - i <- (x$states != noplot) + i <- !(x$states %in% noplot) if (all(i) || !any(i)) { # the !any is a failsafe, in case none are kept then ignore noplot ssurv <- smat(x$pstate) diff -Nru survival-3.1-8/R/print.coxph.R survival-3.1-11/R/print.coxph.R --- survival-3.1-8/R/print.coxph.R 2019-10-28 14:11:43.000000000 +0000 +++ survival-3.1-11/R/print.coxph.R 2020-03-03 19:36:07.000000000 +0000 @@ -36,14 +36,24 @@ # print it group by group tmap <- x$cmap[-1,,drop=FALSE] # ignore the intercept (strata) cname <- colnames(tmap) + printed <- rep(FALSE, length(cname)) for (i in 1:length(cname)) { - tmp2 <- tmp[tmap[,i],, drop=FALSE] - names(dimnames(tmp2)) <- c(cname[i], "") - # restore character row names - rownames(tmp2) <- rownames(tmap)[tmap[,i]>0] - printCoefmat(tmp2, digits=digits, P.values=TRUE, has.Pvalue=TRUE, - signif.stars = signif.stars, ...) - cat("\n") + # if multiple colums of tmat are identical, only print that + # set of coefficients once + if (!printed[i]) { + j <- apply(tmap[-1,, drop=FALSE], 2, + function(x) all(x == tmap[-1,i])) + printed[j] <- TRUE + + tmp2 <- tmp[tmap[,i],, drop=FALSE] + names(dimnames(tmp2)) <- c(paste(cname[j], collapse=", "), "") + # restore character row names + rownames(tmp2) <- rownames(tmap)[tmap[,i]>0] + printCoefmat(tmp2, digits=digits, P.values=TRUE, + has.Pvalue=TRUE, + signif.stars = signif.stars, ...) + cat("\n") + } } cat(" States: ", paste(paste(seq(along=x$states), x$states, sep='= '), diff -Nru survival-3.1-8/R/print.summary.coxph.R survival-3.1-11/R/print.summary.coxph.R --- survival-3.1-8/R/print.summary.coxph.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/R/print.summary.coxph.R 2020-03-05 13:29:54.000000000 +0000 @@ -25,16 +25,50 @@ return() } + if (!is.null(x$cmap)) { # this was a coxphms object + signif.stars <- FALSE #work around issue with printCoefmat + # print it group by group + tmap <- x$cmap[-1,,drop=FALSE] # ignore the intercept (strata) + cname <- colnames(tmap) + printed <- rep(FALSE, length(cname)) + for (i in 1:length(cname)) { + # if multiple colums of tmat are identical, only print that + # set of coefficients once + if (!printed[i]) { + j <- apply(tmap[-1,, drop=FALSE], 2, + function(x) all(x == tmap[-1,i])) + printed[j] <- TRUE + + tmp2 <- x$coefficients[tmap[,i],, drop=FALSE] + names(dimnames(tmp2)) <- c(paste(cname[j], collapse=", "), "") + # restore character row names + rownames(tmp2) <- rownames(tmap)[tmap[,i]>0] + + printCoefmat(tmp2, digits=digits, P.values=TRUE, + has.Pvalue=TRUE, signif.legend=FALSE, + signif.stars = signif.stars, ...) - if(!is.null(x$coefficients)) { - cat("\n") - printCoefmat(x$coefficients, digits=digits, - signif.stars=signif.stars, ...) - } - if(!is.null(x$conf.int)) { - cat("\n") - print(x$conf.int) + if (!is.null(x$conf.int)) { + tmp2 <- x$conf.int[tmap[,i],, drop=FALSE] + rownames(tmp2) <- rownames(tmap)[tmap[,i] >0] + names(dimnames(tmp2)) <- c(paste(cname[j], collapse=", "),"") + print(tmp2, digits=digits, ...) + } + } + } + cat("\n States:", paste(paste(seq(along=x$states), x$states, sep='= '), + collapse=", "), '\n') + } else { + if(!is.null(x$coefficients)) { + cat("\n") + printCoefmat(x$coefficients, digits=digits, + signif.stars=signif.stars, ...) } + if(!is.null(x$conf.int)) { + cat("\n") + print(x$conf.int) + } + } cat("\n") if (!is.null(x$concordance)) { diff -Nru survival-3.1-8/R/print.survcheck.R survival-3.1-11/R/print.survcheck.R --- survival-3.1-8/R/print.survcheck.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/R/print.survcheck.R 2020-03-05 13:29:54.000000000 +0000 @@ -4,18 +4,17 @@ dput(cl) cat("\n") } - #browser() + temp <- x$n + names(temp) <- c("Unique identifiers", "Observations", "Transitions") + print(temp) if(!is.null(x$na.action)){ cat(length(x$na.action), "observations with missing values:","\n") - } - cat(sum(x$events[1,]),"subjects available for analysis","\n") + } - ## change between states - if(nrow(x$transitions)>1){ - cat("Transitions table:\n") - print(x$transitions) - cat('\n') - } + cat("\nTransitions table:\n") + print(x$transitions) + cat('\n') + ## how many of each state does each id have? cat("Number of subjects with 0, 1, ... transitions to each state:\n") print(x$events) @@ -59,70 +58,67 @@ cat("\n") } + temp <- object$n + names(temp) <- c("Unique identifiers", "Observations", "Transitions") + print(temp) + if(!is.null(object$na.action)){ cat(length(object$na.action), "observations with missing values:","\n") } - cat(sum(object$statecount[,1]),"subjects available for analysis","\n") ## change between states - if(nrow(object$transitions)>1){ - cat("Transitions table:\n") - print(object$transitions) - cat('\n') - } + cat("Transitions table:\n") + print(object$transitions) + cat('\n') + ## how many of each state does each id have? cat("Number of subjects with 1, 2, ... copies of each state:\n") print(object$events) cat("\n") - ## change between states - if(nrow(object$transitions)>1){ - cat("# subjects moving between states:\n") - print(object$transitions) - cat('\n') - } if(object$flag["overlap"]>0) { - cat("Overlap check: ", + cat("Overlap: ", length(object$overlap$id), ifelse(length(object$overlap$id)==1," id (", " ids ("), length(object$overlap$row), " rows)\n", sep="") - o.mat <- cbind(id=object$id,object$Y)[object$overlap$row,] - max.show <- min(max.show, length(object$overlap$id)) - print(o.mat[o.mat[,1] %in% object$overlap$id[1:max.show],]) + odat <- data.frame(id=object$id, y= object$Y)[object$overlap$row,] + nshow <- min(max.show, length(object$overlap$id)) + print(odat[odat[,1] %in% object$overlap$id[1:nshow],]) cat('\n') } if(object$flag["gap"]>0) { - cat("Gap check: ", + cat("Gap: ", length(object$gap$id), ifelse(length(object$gap$id)==1," id (", " ids ("), length(object$gap$row), " rows)\n", sep="") - g.mat <- cbind(id=object$id,object$Y)[object$gap$row,] - max.show <- min(max.show, length(object$gap$id)) - print(g.mat[g.mat[,1] %in% object$gap$id[1:max.show],]) + gdat <- data.frame(id=object$id, y=object$Y)[object$gap$row,] + nshow <- min(max.show, length(object$gap$id)) + print(gdat[gdat[,1] %in% object$gap$id[1:nshow],]) cat('\n') } + if(object$flag["teleport"]>0) { - cat("Teleport check: ", + cat("Teleport: ", length(object$teleport$id), ifelse(length(object$teleport$id)==1," id (", " ids ("), length(object$teleport$row), " rows)\n", sep="") - t.mat <- cbind(id=object$id,object$Y)[object$teleport$row,] - max.show <- min(max.show, length(object$teleport$id)) - print(t.mat[t.mat[,1] %in% object$teleport$id[1:max.show],]) + tdat <- data.frame(id=object$id,object$Y)[object$teleport$row,] + nshow <- min(max.show, length(object$teleport$id)) + print(tdat[tdat[,1] %in% object$teleport$id[1:nshow],]) cat('\n') } if(object$flag["jump"]){ - cat("Jump check: ", + cat("Jump: ", length(object$jump$id), ifelse(length(object$jump$id)==1," id (", " ids ("), length(object$jump$row), " rows)\n", sep="") - j.mat <- cbind(id=object$id,object$Y)[object$jump$row,] - max.show <- min(max.show, length(object$jump$id)) - print(j.mat[j.mat[,1] %in% object$jump$id[1:max.show],]) + jdat <- data.frame(id=object$id, y= object$Y)[object$jump$row,] + nshow <- min(max.show, length(object$jump$id)) + print(jdat[jdat[,1] %in% object$jump$id[1:nshow],]) cat('\n') } } diff -Nru survival-3.1-8/R/pspline.R survival-3.1-11/R/pspline.R --- survival-3.1-8/R/pspline.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/R/pspline.R 2020-02-14 22:50:48.000000000 +0000 @@ -74,12 +74,14 @@ if (any (combine != floor(combine) | combine < 0) || any(diff(combine) < 0)) stop("combine must be an increasing vector of positive integers") - if (!intercept) combine <- c(0, combine) - if (length(combine) != ncol(newx)) + if (!intercept) ctemp <- c(0, combine) + else ctemp <- combine + # the intercept is removed from newx later + if (length(ctemp) != ncol(newx)) stop("wrong length for combine") - uc <- sort(unique(combine)) + uc <- sort(unique(ctemp)) tmat <- matrix(0., nrow=ncol(newx), ncol=length(uc)) - for (i in 1:length(uc)) tmat[combine==uc[i], i] <- 1 + for (i in 1:length(uc)) tmat[ctemp==uc[i], i] <- 1 newx <- newx %*% tmat } @@ -99,6 +101,7 @@ attributes(newx) <- c(attributes(newx), list(intercept=intercept, nterm=nterm, Boundary.knots=Boundary.knots)) + if (!missing(combine)) attr(newx, "combine") <- combine class(newx) <- "pspline" return(newx) } @@ -188,6 +191,8 @@ attributes(newx) <- c(attributes(newx), temp, list(intercept=intercept, nterm=nterm, Boundary.knots=Boundary.knots)) + if (!missing(combine)) attr(newx, "combine") <- combine + class(newx) <- c("pspline", 'coxph.penalty') newx } @@ -195,16 +200,21 @@ makepredictcall.pspline <- function(var, call) { if (call[[1]] != as.name("pspline")) return(call) #wrong phone number newcall <- call[1:2] #don't let the user override anything - at <- attributes(var)[c("nterm", "intercept", "Boundary.knots", "combine")] + indx <- match(c("nterm", "intercept", "Boundary.knots", "combine"), + names(attributes(var)), nomatch=0) + at <- attributes(var)[indx] newcall[names(at)] <- at + newcall } predict.pspline <- function(object, newx, ...) { if (missing(newx)) return(object) - a <- c(list(x=newx, penalty=FALSE), - attributes(object)[c("intercept, Boundary.knots", "combine")]) - do.call("pspline", a) + indx <- match(c("nterm", "intercept", "Boundary.knots", "combine"), + names(attributes(object)), nomatch=0) + at <- c(list(x=newx, penalty=FALSE), + attributes(object)[indx]) + do.call("pspline", at) } # Given a pspline basis, recover x diff -Nru survival-3.1-8/R/residuals.survreg.R survival-3.1-11/R/residuals.survreg.R --- survival-3.1-8/R/residuals.survreg.R 2019-12-01 17:13:45.000000000 +0000 +++ survival-3.1-11/R/residuals.survreg.R 2020-03-05 13:32:42.000000000 +0000 @@ -155,7 +155,9 @@ else score <- cbind(score, deriv[,4]) } rr <- score %*% vv - if (type=='dfbetas') rr <- rr %*% diag(1/sqrt(diag(vv))) + # cause column names to be retained + # old: if (type=='dfbetas') rr[] <- rr %*% diag(1/sqrt(diag(vv))) + if (type=='dfbetas') rr <- rr * rep(1/sqrt(diag(vv)), each=nrow(rr)) if (type=='ldcase') rr<- rowSums(rr*score) } diff -Nru survival-3.1-8/R/summary.coxph.R survival-3.1-11/R/summary.coxph.R --- survival-3.1-8/R/summary.coxph.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/R/summary.coxph.R 2020-03-05 13:29:54.000000000 +0000 @@ -64,7 +64,11 @@ rval$concordance <- cox$concordance[6:7] names(rval$concordance) <- c("C", "se(C)") } - + if (inherits(cox, "coxphms")) { + rval$cmap <- cox$cmap + rval$states <- cox$states + } + class(rval) <-"summary.coxph" rval } diff -Nru survival-3.1-8/R/survcheck.R survival-3.1-11/R/survcheck.R --- survival-3.1-8/R/survcheck.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/R/survcheck.R 2020-03-05 13:29:54.000000000 +0000 @@ -29,7 +29,8 @@ if (!is.null(istate) && length(istate) !=n) stop("wrong length for istate") fit <- survcheck2(Y, id, istate, istate0) -# fit$istate <- NULL # used by coxph, but not part of user printout + fit$n <- c(id = length(unique(id)), observations =length(id), + transitions = sum(fit$transitions)) na.action <- attr(mf, "na.action") if (!is.null(na.action)) { fit$na.action <- na.action @@ -44,7 +45,8 @@ } } } - fit$Y <- Y # used by the summary function + + fit$Y <- Y # used by the summary function for details fit$id <- unname(id) # "" "" fit$call <- Call class(fit) <- "survcheck" @@ -67,8 +69,8 @@ survcheck2 <- function(y, id, istate=NULL, istate0="(s0)") { n <- length(id) ny <- ncol(y) - # the next few line are a debug for my code, since survcheck2 is not visible - # to users + # the next few line are a debug for my code; survcheck2 is not visible + # to users so only survival can call it directly if (!is.Surv(y) || is.null(attr(y, "states")) || any(y[,ncol(y)] > length(attr(y, "states")))) stop("survcheck2 called with an invalid y argument") @@ -80,7 +82,8 @@ } else { if (length(istate) !=n) stop ("wrong length for istate") - cstate <- as.factor(istate) + if (is.factor(istate)) cstate <- istate[, drop=TRUE] #drop unused levels + else cstate <- as.factor(istate) inull <- FALSE } @@ -92,55 +95,69 @@ index <- match(levels(cstate), ystate, nomatch=0) states <- c(levels(cstate)[index==0], ystate) cstate2 <- factor(cstate, states) - # we keep a form with all the levels for returning to the parent (cstate2) - # one without this (cstate) to make a smaller transitions table # Calculate the counts per id for each state, e.g., 10 subjects had # 3 visits to state 2, etc. - # Don't count "censored" as an endpoint, nor any missings. But we can't - # just omit censored rows, or those with 0 transitions don't get counted! - # Instead remove the 'censored' column after making tab1 + # Count the censors, so that each subject gets a row in the table, + # but then toss that column tab1 <- table(id, factor(y[,ncol(y)], 0:length(ystate)))[,-1, drop=FALSE] tab1 <- cbind(tab1, rowSums(tab1)) tab1.levels <- sort(unique(c(tab1))) #unique counts events <- apply(tab1, 2, function(x) table(factor(x, tab1.levels))) dimnames(events) = list("count"= tab1.levels, "state"= c(ystate, "(any)")) - - # check for errors + + # Use a C routine to create 3 variables: a: an index of whether this is + # the first (1) or last(2) observation for a subject, 3=both, 0=neither, + # b. current state, and + # c. sign of (start of this interval - end of prior one) + # start by making stat2 = status re-indexed to the full set of states + ny <- ncol(y) sindx <- match(ystate, states) - if (ncol(y)==2) y <- cbind(0,y) # make it 3 cols for the C routine - stat2 <- ifelse(y[,3]==0, 0L, - sindx[pmax(1, y[,3])]) # map the status - index <- order(id, y[,ny-1], y[,1]) # by start time, within stop time - id2 <- match(id, id) # we need unique integers, but don't need 1, 2,... - check <- .Call(Cmulticheck, y, stat2, id2, as.integer(cstate2), index-1L) - if (inull) { + stat2 <- ifelse(y[,ny]==0, 0L, sindx[pmax(1L, y[,ny])]) + id2 <- match(id, unique(id)) # we need unique integers + if (ncol(y)==2) { + index <- order(id, y[,1]) + check <- .Call(Cmulticheck, rep(0., n), y[,1], stat2, id2, + as.integer(cstate2), index- 1L) + } else { + index <- order(id, y[,2], y[,1]) + check <- .Call(Cmulticheck, y[,1], y[,2], stat2, id2, + as.integer(cstate2), index- 1L) + } + + if (inull && ny> 2) { # if there was no istate entered in, use the constructed one from # the check routine + # if ny=2 then every row starts at time 0 cstate2 <-factor(check$cstate, seq(along=states), states) } # create the transtions table # if someone has an intermediate visit, i.e., (0,10, 0)(10,20,1), don't # report the false 'censoring' in the transitions table - yfac <- factor(y[,3], c(seq(along=ystate), 0), to.names) - keep <- (y[index,3]!=0 | !duplicated(id[index], fromLast=TRUE)) - keep[index] <- keep # put it back into data order + # make it compact by removing any rows or cols that are all 0 + keep <- (stat2 !=0 | check$dupid > 1) # not censored or last obs of this id transitions <- table(from=cstate2[keep], - to= yfac[keep], + to= factor(stat2[keep], c(seq(along=states), 0), + c(states, "(censored)")), useNA="ifany") + transitions <- transitions[rowSums(transitions) >0, colSums(transitions)>0, + drop = FALSE] - # now continue with error checks. + # now continue with error checks + # A censoring hole in the middle, such as happens with survSplit, + # uses "last state carried forward" in Cmultistate, which also + # sets the "gap" to 0 for the first obs of a subject mismatch <- (as.numeric(cstate2) != check$cstate) - + # gap = 0 (0, 10], (10, 15] # gap = 1 (0, 10], (12, 15] # a hole in the time # gap = -1 (0, 10], (9, 15] # overlapping times flag <- c(overlap= sum(check$gap < 0), gap = sum(check$gap > 0 & !mismatch), jump = sum(check$gap > 0 & mismatch), - teleport = sum(check$gap==0 & mismatch & check$dupid==1)) + teleport = sum(check$gap==0 & mismatch & check$dupid%%2 ==0)) rval <- list(states=states, transitions=transitions, events= t(events), flag=flag, @@ -160,7 +177,7 @@ rval$jump <- list(row=j, id= unique(id[j])) } if (flag["teleport"] > 0) { - j <- (check$gap==0 & mismatch) + j <- which(check$gap==0 & mismatch) rval$teleport <- list(row=j, id= unique(id[j])) } diff -Nru survival-3.1-8/R/tmerge.R survival-3.1-11/R/tmerge.R --- survival-3.1-8/R/tmerge.R 2019-12-01 17:13:50.000000000 +0000 +++ survival-3.1-11/R/tmerge.R 2020-03-05 13:32:47.000000000 +0000 @@ -317,8 +317,8 @@ if (argclass[ii] == "tdc") { if (argname[[ii]] %in% tevent) stop("attempt to turn event variable", argname[[ii]], "into a tdc") - #if (!is.null(newvar)) - #warning(paste0("replacement of variable '", argname[ii], "'")) + if (!is.null(newvar)) + warning(paste0("replacement of variable '", argname[ii], "'")) # id can be any data type; feed integers to the C routine storage.mode(dstop) <- storage.mode(etime) <- "double" #if time is integer diff -Nru survival-3.1-8/src/init.c survival-3.1-11/src/init.c --- survival-3.1-8/src/init.c 2019-09-30 11:53:47.000000000 +0000 +++ survival-3.1-11/src/init.c 2020-03-05 13:29:54.000000000 +0000 @@ -52,7 +52,7 @@ {"Cgchol", (DL_FUNC) &gchol, 2}, {"Cgchol_solve", (DL_FUNC) &gchol_solve, 3}, {"Cgchol_inv", (DL_FUNC) &gchol_inv, 2}, - {"Cmulticheck", (DL_FUNC) &multicheck, 5}, + {"Cmulticheck", (DL_FUNC) &multicheck, 6}, {"Cpyears3b", (DL_FUNC) &pyears3b, 10}, {"Csurvfitci", (DL_FUNC) &survfitci, 11}, {"Csurvfitkm", (DL_FUNC) &survfitkm, 9}, diff -Nru survival-3.1-8/src/multicheck.c survival-3.1-11/src/multicheck.c --- survival-3.1-8/src/multicheck.c 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/src/multicheck.c 2020-03-05 15:22:55.000000000 +0000 @@ -1,17 +1,19 @@ /* ** create three variables -** 1. dupid = 1 if the id is the second or more +** 1. dupid = 1 first obs for this id, 2=last, 3=both, 0=neither ** 2. gap = -1 the start time for this obs is before the prior end time ** 1 the start time for this obs is after the prior end time ** 0 times match, or the first instance of the observation ** 3. cstate = current state for this obs, obtained by chaining forward -** initial state at entry, first transtion, second, etc. +** initial state at entry, first transtion, second, etc. but ignoring +** censors as a 'change'. */ #include "survS.h" #include "survproto.h" -SEXP multicheck(SEXP y2, SEXP status2, SEXP id2, SEXP istate2, SEXP sort2) { - int i, ii, ilag; +SEXP multicheck(SEXP time12, SEXP time22, SEXP status2, SEXP id2, + SEXP istate2, SEXP sort2) { + int i, ii, oldid, oldii; int n; double *time1, *time2; @@ -25,8 +27,8 @@ SEXP rlist; /* return list*/ n = LENGTH(id2); - time1 = REAL(y2); - time2 = time1 + n; + time1 = REAL(time12); + time2 = REAL(time22); status = INTEGER(status2); id = INTEGER(id2); istate = INTEGER(istate2); @@ -37,27 +39,31 @@ gap = INTEGER(SET_VECTOR_ELT(rlist, 1, allocVector(INTSXP, n))); cstate= INTEGER(SET_VECTOR_ELT(rlist, 2, allocVector(INTSXP, n))); - ilag = -1; + oldid = -1; /* the input id values are all > 0 */ + oldii = 0 ; /* stop a complaint from -Wall option of gcc */ cstate[0] = istate[0]; for (i=0; i time2[ilag]) gap[ii] =1; + if (id[ii] == oldid) { + dupid[ii] = 0; + if (time1[ii] == time2[oldii]) gap[ii] =0; + else if (time1[ii] > time2[oldii]) gap[ii] =1; else gap[ii] = -1; - if (status[ilag]>0) cstate[ii] = status[ilag]; - else cstate[ii] = cstate[ilag]; + if (status[oldii]>0) cstate[ii] = status[oldii]; + else cstate[ii] = cstate[oldii]; } else { /* first obs for a new id*/ + oldid = id[ii]; dupid[ii] =0; gap[ii] =0; cstate[ii] = istate[ii]; + if (i>0) dupid[oldii] += 2; /* prior obs was last for that id*/ } - ilag = ii; + oldii = ii; } + dupid[oldii] += 2; UNPROTECT(1); return(rlist); } diff -Nru survival-3.1-8/src/survproto.h survival-3.1-11/src/survproto.h --- survival-3.1-8/src/survproto.h 2019-12-01 16:48:18.000000000 +0000 +++ survival-3.1-11/src/survproto.h 2020-03-05 13:29:54.000000000 +0000 @@ -156,7 +156,8 @@ void init_doloop(int min, int max); int doloop (int nloops, int *index); -SEXP multicheck(SEXP y2, SEXP status2, SEXP id2, SEXP istate2, SEXP sort2); +SEXP multicheck(SEXP time12, SEXP time22, SEXP status2, SEXP id2, + SEXP istate2, SEXP sort2); int *norisk(int n, double *time1, double *time2, double *status, int *sort1, int *sort2, int *strata); diff -Nru survival-3.1-8/src/zph1.c survival-3.1-11/src/zph1.c --- survival-3.1-8/src/zph1.c 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/src/zph1.c 2020-02-28 12:14:50.000000000 +0000 @@ -31,6 +31,7 @@ double *a, *a2, **cmat, **cmat2; double denom=0, dtime, ndead, denom2; double risk, deadwt, wtave; + int nprotect; /* scalar input arguments and counts*/ int nused, nvar, nevent, nstrat; @@ -74,10 +75,14 @@ /* ** Set up the ragged arrays and scratch space */ - if (MAYBE_REFERENCED(covar2)) - covar = dmatrix(REAL(duplicate(covar2)), nused, nvar); + nprotect =0; + if (MAYBE_REFERENCED(covar2)){ + nprotect++; + covar = dmatrix(REAL(PROTECT(duplicate(covar2))), nused, nvar); + } else covar = dmatrix(REAL(covar2), nused, nvar); + nprotect++; PROTECT(rlist = mkNamed(VECSXP, rnames)); u = REAL(SET_VECTOR_ELT(rlist, 0, allocVector(REALSXP, 2*nvar))); dtemp = REAL(SET_VECTOR_ELT(rlist, 1, allocMatrix(REALSXP, 2*nvar, 2*nvar))); @@ -294,6 +299,6 @@ imat[i+nvar][j] = imat[j][i+nvar]; } - UNPROTECT(1); + UNPROTECT(nprotect); return(rlist); } diff -Nru survival-3.1-8/src/zph2.c survival-3.1-11/src/zph2.c --- survival-3.1-8/src/zph2.c 2019-10-31 12:04:55.000000000 +0000 +++ survival-3.1-11/src/zph2.c 2020-02-28 12:13:41.000000000 +0000 @@ -19,6 +19,7 @@ int *keep; double denom, dtime=0, ndead, denom2; double risk, meanwt; + int nprotect; /* scalar input arguments and counts*/ int nused, nvar, nevent, nstrat; @@ -63,10 +64,14 @@ /* ** Set up the ragged arrays and scratch space */ - if (MAYBE_REFERENCED(covar2)) - covar = dmatrix(REAL(duplicate(covar2)), nused, nvar); + nprotect =0; + if (MAYBE_REFERENCED(covar2)) { + nprotect =1; + covar = dmatrix(REAL(PROTECT(duplicate(covar2))), nused, nvar); + } else covar = dmatrix(REAL(covar2), nused, nvar); + nprotect++; PROTECT(rlist = mkNamed(VECSXP, rnames)); u = REAL(SET_VECTOR_ELT(rlist, 0, allocVector(REALSXP, 2*nvar))); dtemp = REAL(SET_VECTOR_ELT(rlist, 1, allocMatrix(REALSXP, 2*nvar, 2*nvar))); @@ -363,6 +368,6 @@ imat[i+nvar][j] = imat[j][i+nvar]; } - UNPROTECT(1); + UNPROTECT(nprotect); return(rlist); } diff -Nru survival-3.1-8/tests/book5.Rout.save survival-3.1-11/tests/book5.Rout.save --- survival-3.1-8/tests/book5.Rout.save 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/tests/book5.Rout.save 2020-02-28 13:29:24.000000000 +0000 @@ -1,6 +1,6 @@ -R Under development (unstable) (2019-01-23 r76006) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing +R Under development (unstable) (2020-02-28 r77869) -- "Unsuffered Consequences" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -162,7 +162,7 @@ exp(coef) exp(-coef) lower .95 upper .95 x 2.362 0.4233 0.5839 9.556 -Concordance= 0.638 (se = 0.161 ) +Concordance= 0.637 (se = 0.161 ) Likelihood ratio test= 1.69 on 1 df, p=0.2 Wald test = 1.45 on 1 df, p=0.2 Score (logrank) test = 1.52 on 1 df, p=0.2 @@ -203,4 +203,4 @@ > > proc.time() user system elapsed - 0.860 0.040 0.902 + 0.906 0.043 0.941 diff -Nru survival-3.1-8/tests/book6.Rout.save survival-3.1-11/tests/book6.Rout.save --- survival-3.1-8/tests/book6.Rout.save 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/tests/book6.Rout.save 2020-03-05 15:17:18.000000000 +0000 @@ -1,6 +1,6 @@ -R Under development (unstable) (2019-01-23 r76006) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing +R Under development (unstable) (2020-02-28 r77869) -- "Unsuffered Consequences" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -151,7 +151,7 @@ exp(coef) exp(-coef) lower .95 upper .95 x 2.393 0.4179 0.5921 9.672 -Concordance= 0.638 (se = 0.161 ) +Concordance= 0.637 (se = 0.161 ) Likelihood ratio test= 1.75 on 1 df, p=0.2 Wald test = 1.5 on 1 df, p=0.2 Score (logrank) test = 1.58 on 1 df, p=0.2 @@ -187,4 +187,4 @@ > > proc.time() user system elapsed - 0.832 0.032 0.875 + 0.891 0.068 0.952 diff -Nru survival-3.1-8/tests/doweight.Rout.save survival-3.1-11/tests/doweight.Rout.save --- survival-3.1-8/tests/doweight.Rout.save 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/tests/doweight.Rout.save 2020-02-28 13:29:52.000000000 +0000 @@ -1,6 +1,6 @@ -R Under development (unstable) (2019-05-16 r76518) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing +R Under development (unstable) (2020-02-28 r77869) -- "Unsuffered Consequences" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -89,7 +89,7 @@ exp(coef) exp(-coef) lower .95 upper .95 x 2.362 0.4233 0.5839 9.556 -Concordance= 0.638 (se = 0.161 ) +Concordance= 0.637 (se = 0.161 ) Likelihood ratio test= 1.69 on 1 df, p=0.2 Wald test = 1.45 on 1 df, p=0.2 Score (logrank) test = 1.52 on 1 df, p=0.2 @@ -200,7 +200,7 @@ exp(coef) exp(-coef) lower .95 upper .95 x 2.362 0.4233 0.5839 9.556 -Concordance= 0.638 (se = 0.172 ) +Concordance= 0.637 (se = 0.172 ) Likelihood ratio test= 1.69 on 1 df, p=0.2 Wald test = 1.45 on 1 df, p=0.2 Score (logrank) test = 1.52 on 1 df, p=0.2 @@ -422,4 +422,4 @@ > > proc.time() user system elapsed - 0.829 0.065 0.888 + 0.986 0.049 1.028 diff -Nru survival-3.1-8/tests/r_sas.Rout.save survival-3.1-11/tests/r_sas.Rout.save --- survival-3.1-8/tests/r_sas.Rout.save 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/tests/r_sas.Rout.save 2020-02-28 12:35:35.000000000 +0000 @@ -1,6 +1,6 @@ -R Under development (unstable) (2019-03-27 r76274) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing +R Under development (unstable) (2020-01-13 r77659) -- "Unsuffered Consequences" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -424,10 +424,10 @@ > points(kfit$time, qnorm(1-kfit$surv), pch='+') > > dev.off() #close the plot file -pdf - 2 +null device + 1 > > > proc.time() user system elapsed - 2.414 0.067 2.476 + 2.655 0.063 2.712 diff -Nru survival-3.1-8/tests/sexpm.save survival-3.1-11/tests/sexpm.save --- survival-3.1-8/tests/sexpm.save 1970-01-01 00:00:00.000000000 +0000 +++ survival-3.1-11/tests/sexpm.save 2020-03-05 13:29:54.000000000 +0000 @@ -0,0 +1,59 @@ +# Test the sexpm functions that are used for fast matrix exponentials +# +library(Matrix) +library(survival) +aeq <- function(x, y, ...) all.equal(as.vector(x),as.vector(y), ...) + +nfun <- length(sexpm) +times <- c(.1, 1, 5) +test1 <- matrix(FALSE, nfun, length(times), + dimnames=list(names(sexpm), paste0("t=", times))) +test2 <- test1 +eps <- 1e-8 + +dtest <- function(x, tmat, time= 1.3, eps=1e-8) { + # Check a derivative + if (missing(tmat)) { + tmat <- matrix(0, x$nstate, x$nstate) + tmat[x$nonzero] <- runif(length(x$nonzero), .1, 3) + diag(tmat) <- -rowSums(tmat) + } + else { + if (!aeq(dim(tmat), c(x$nstate, x$nstate))) stop("invalid tmat") + diag(tmat) <- diag(tmat) - rowSums(tmat) + } + + d1 <- x$deriv(tmat, time) + d2 <- 0*d1 + nz <- x$nonzero + exp1 <- as.matrix(expm(tmat* time)) + + for (j in 1:length(nz)) { + temp <- tmat + temp[nz[j]] <- temp[nz[j]] + eps + diag(temp) <- diag(temp) - rowSums(temp) + exp2 <- as.matrix(expm(temp* time)) + d2[,,j] <- (exp2 - exp1)/eps + } + list(d1=d1, d2=d2) +} + +for (i in 1:nfun) { + j <- sexpm[[i]] + tmat <- matrix(0, j$nstate, j$nstate) + tmat[j$nonzero] <- runif(length(j$nonzero), .1, 4) + diag(tmat) <- -rowSums(tmat) + + dtemp <- logical(length(j$nonzero)) + for (k in 1:length(times)) { + m1 <- j$mexp(tmat, times[k]) + m2 <- expm(tmat*times[k]) + test1[i,k] <- isTRUE(aeq(m1, as.matrix(m2))) + # now derivatives + temp <- dtest(j, tmat, times[k], eps) + test2[i,k] <- isTRUE(aeq(temp$d1, temp$d2, tol=sqrt(eps))) + } +} + +all(test1) # should all be TRUE +all(test2) diff -Nru survival-3.1-8/tests/survcheck.R survival-3.1-11/tests/survcheck.R --- survival-3.1-8/tests/survcheck.R 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/tests/survcheck.R 2020-03-05 15:39:55.000000000 +0000 @@ -20,7 +20,8 @@ t1 = c(0, 10, 20, 0, 15,25, 10, 15, 18), t2 = c(10, 20,30, 20, 24, 26, 13, 18, 25), status= factor(c(0, 1, 2, 1,2, 0, 1, 0, 3)), - x = c(1:5, NA, 7:9)) + x = c(1:5, NA, 7:9), + stringsAsFactors = FALSE) fit2 <- survcheck(Surv(t1, t2, status) ~ 1, data2, id=id) aeq(fit2$flag , c(1,2,0,0)) @@ -38,24 +39,26 @@ # let a missing value in fit2b <- survcheck(Surv(t1, t2, status) ~ x, data2, id=id) aeq(fit2b$flag , c(1,1,0,0)) -aeq(fit2b$transition, rbind(c(3,0,0,0), c(0,2,1,0))) +aeq(fit2b$transition, rbind(c(3,0,0), c(0,2,1))) (fit2b$overlap$id == 'B') (fit2b$overlap$row ==5) all(fit2b$gap$id == "C") aeq(fit2b$gap$row, 8) +# designed to trigger all 4 error types data3 <- data2 levels(data3$status) <- c("cens", "mgus", "recur", "death") data3$istate <- c("entry", "entry", "recur", "entry", "recur", "recur", "entry", "recur", "recur") fit3 <- survcheck(Surv(t1, t2, status) ~ 1, data3, id=id, istate=istate) -aeq(fit3$flag, c(1, 0, 2, 1)) +aeq(fit3$flag, c(1, 1, 1, 2)) aeq(fit3$transitions, rbind(c(3,0,0,0), c(0,2,1,1))) all.equal(fit3$overlap, fit2$overlap) all(fit3$teleport$id == c("A", "C")) all(fit3$teleport$row == c(3,9)) all(fit3$jump$id == "C") all(fit3$jump$row == 8) +all.equal(fit3$gap, list(row=6L, id= "B")) diff -Nru survival-3.1-8/tests/survcheck.Rout.save survival-3.1-11/tests/survcheck.Rout.save --- survival-3.1-8/tests/survcheck.Rout.save 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/tests/survcheck.Rout.save 2020-03-05 15:45:34.000000000 +0000 @@ -1,6 +1,6 @@ -R Under development (unstable) (2019-08-23 r77061) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing +R Under development (unstable) (2020-02-28 r77869) -- "Unsuffered Consequences" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -29,7 +29,7 @@ > aeq(fit1$overlap$id, 2) [1] TRUE > aeq(fit1$transitions, c(2,1,1,1)) -[1] "Numeric: lengths (16, 4) differ" +[1] TRUE > > > # dummy data 2, no initial values, start stop data @@ -41,13 +41,14 @@ + t1 = c(0, 10, 20, 0, 15,25, 10, 15, 18), + t2 = c(10, 20,30, 20, 24, 26, 13, 18, 25), + status= factor(c(0, 1, 2, 1,2, 0, 1, 0, 3)), -+ x = c(1:5, NA, 7:9)) ++ x = c(1:5, NA, 7:9), ++ stringsAsFactors = FALSE) > fit2 <- survcheck(Surv(t1, t2, status) ~ 1, data2, id=id) > > aeq(fit2$flag , c(1,2,0,0)) [1] TRUE > aeq(fit2$transition, rbind(c(3,0,0,0), c(0,2,1,0), c(0,0,0,1))) -[1] "Numeric: lengths (16, 12) differ" +[1] TRUE > (fit2$overlap$id == 'B') [1] TRUE > (fit2$overlap$row ==5) @@ -67,8 +68,8 @@ > fit2b <- survcheck(Surv(t1, t2, status) ~ x, data2, id=id) > aeq(fit2b$flag , c(1,1,0,0)) [1] TRUE -> aeq(fit2b$transition, rbind(c(3,0,0,0), c(0,2,1,0))) -[1] "Numeric: lengths (16, 8) differ" +> aeq(fit2b$transition, rbind(c(3,0,0), c(0,2,1))) +[1] TRUE > (fit2b$overlap$id == 'B') [1] TRUE > (fit2b$overlap$row ==5) @@ -78,32 +79,32 @@ > aeq(fit2b$gap$row, 8) [1] TRUE > +> # designed to trigger all 4 error types > data3 <- data2 > levels(data3$status) <- c("cens", "mgus", "recur", "death") > data3$istate <- c("entry", "entry", "recur", "entry", "recur", "recur", + "entry", "recur", "recur") > fit3 <- survcheck(Surv(t1, t2, status) ~ 1, data3, id=id, istate=istate) > -> aeq(fit3$flag, c(1, 0, 2, 1)) -[1] "Mean relative difference: 0.75" +> aeq(fit3$flag, c(1, 1, 1, 2)) +[1] TRUE > aeq(fit3$transitions, rbind(c(3,0,0,0), c(0,2,1,1))) -[1] "Numeric: lengths (16, 8) differ" +[1] TRUE > all.equal(fit3$overlap, fit2$overlap) [1] TRUE > all(fit3$teleport$id == c("A", "C")) [1] TRUE > all(fit3$teleport$row == c(3,9)) -[1] FALSE -Warning message: -In fit3$teleport$row == c(3, 9) : - longer object length is not a multiple of shorter object length +[1] TRUE > all(fit3$jump$id == "C") [1] TRUE > all(fit3$jump$row == 8) [1] TRUE +> all.equal(fit3$gap, list(row=6L, id= "B")) +[1] TRUE > > > > proc.time() user system elapsed - 0.993 0.059 1.163 + 0.817 0.051 0.862 diff -Nru survival-3.1-8/tests/tmerge3.R survival-3.1-11/tests/tmerge3.R --- survival-3.1-8/tests/tmerge3.R 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/tests/tmerge3.R 2020-02-28 12:33:00.000000000 +0000 @@ -6,9 +6,9 @@ # a row merged in that is a censor, don't mark it as a cumevent. # base <- data.frame( - id = 1:2, tstart = c(0, 0), tstop = c(10, 10), got_flue = c(0, 0), - has_flue = factor(c("no", "no"), levels = c("no", "yes"))) -base <- tmerge(base, base, id = id, got_flue = event(tstop, got_flue)) + id = 1:2, tstart = c(0, 0), tstop = c(10, 10), got_flu = c(0, 0), + has_flu = factor(c("no", "no"), levels = c("no", "yes"))) +base <- tmerge(base, base, id = id, got_flu = event(tstop, got_flu)) # add time-varying covariates vars <- data.frame(id = c(1, rep(2, 5)), time = c(0, (0:4) * 2), x = rnorm(6)) @@ -18,13 +18,13 @@ events <- data.frame( id = c(2, 2, 2), # notice the zero -- the second row should not add an event - got_flue = c(1,0,2), - has_flue = c("yes", "no", "yes"), + got_flu = c(1,0,2), + has_flu = c("yes", "no", "yes"), time = c(3, 5, 8)) -b2 <- tmerge(base, events, id = id, got_flue = cumevent(time, got_flue), - has_flue = tdc(time, has_flue)) +b2 <- tmerge(base, events, id = id, got_flu = cumevent(time, got_flu), + has_flu = tdc(time, has_flu)) -all.equal(b2$got_flue, c(0,0,1,0,0,0,3,0)) +all.equal(b2$got_flu, c(0,0,1,0,0,0,3,0)) # Tied times in the merger data set diff -Nru survival-3.1-8/tests/tmerge3.Rout.save survival-3.1-11/tests/tmerge3.Rout.save --- survival-3.1-8/tests/tmerge3.Rout.save 2019-09-16 12:18:11.000000000 +0000 +++ survival-3.1-11/tests/tmerge3.Rout.save 2020-02-28 12:40:20.000000000 +0000 @@ -1,6 +1,6 @@ -R Under development (unstable) (2019-08-23 r77061) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing +R Under development (unstable) (2020-01-13 r77659) -- "Unsuffered Consequences" +Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. @@ -23,9 +23,9 @@ > # a row merged in that is a censor, don't mark it as a cumevent. > # > base <- data.frame( -+ id = 1:2, tstart = c(0, 0), tstop = c(10, 10), got_flue = c(0, 0), -+ has_flue = factor(c("no", "no"), levels = c("no", "yes"))) -> base <- tmerge(base, base, id = id, got_flue = event(tstop, got_flue)) ++ id = 1:2, tstart = c(0, 0), tstop = c(10, 10), got_flu = c(0, 0), ++ has_flu = factor(c("no", "no"), levels = c("no", "yes"))) +> base <- tmerge(base, base, id = id, got_flu = event(tstop, got_flu)) > > # add time-varying covariates > vars <- data.frame(id = c(1, rep(2, 5)), time = c(0, (0:4) * 2), x = rnorm(6)) @@ -35,13 +35,16 @@ > events <- data.frame( + id = c(2, 2, 2), + # notice the zero -- the second row should not add an event -+ got_flue = c(1,0,2), -+ has_flue = c("yes", "no", "yes"), ++ got_flu = c(1,0,2), ++ has_flu = c("yes", "no", "yes"), + time = c(3, 5, 8)) -> b2 <- tmerge(base, events, id = id, got_flue = cumevent(time, got_flue), -+ has_flue = tdc(time, has_flue)) +> b2 <- tmerge(base, events, id = id, got_flu = cumevent(time, got_flu), ++ has_flu = tdc(time, has_flu)) +Warning message: +In tmerge(base, events, id = id, got_flu = cumevent(time, got_flu), : + replacement of variable 'has_flu' > -> all.equal(b2$got_flue, c(0,0,1,0,0,0,3,0)) +> all.equal(b2$got_flu, c(0,0,1,0,0,0,3,0)) [1] TRUE > > @@ -77,4 +80,4 @@ > > proc.time() user system elapsed - 0.702 0.051 0.747 + 0.808 0.035 0.836 diff -Nru survival-3.1-8/vignettes/survival.Rnw survival-3.1-11/vignettes/survival.Rnw --- survival-3.1-8/vignettes/survival.Rnw 2019-10-28 15:08:55.000000000 +0000 +++ survival-3.1-11/vignettes/survival.Rnw 2020-01-28 21:37:52.000000000 +0000 @@ -460,10 +460,10 @@ Other options: \begin{itemize} - \item xaxt('r') It has been traditional to have survival curves touch + \item xaxs('r') It has been traditional to have survival curves touch the left axis (I will not speculate as to why). - This can be accomplished using \code{xaxt='S'}, which was the default in + This can be accomplished using \code{xaxs='S'}, which was the default before survival 3.x. The current default is the standard R style, which leaves space between the curve and the axis. \item The follow-up time in the data set is in days. This is very common in diff -Nru survival-3.1-8/vignettes/timedep.Rnw survival-3.1-11/vignettes/timedep.Rnw --- survival-3.1-8/vignettes/timedep.Rnw 2019-09-16 12:17:55.000000000 +0000 +++ survival-3.1-11/vignettes/timedep.Rnw 2020-02-06 17:04:24.000000000 +0000 @@ -278,7 +278,7 @@ taken from \code{data2}. The idea behind the function is that each addition will be ``slipped in'' to the original data in the same way that one -would slide a new card into an existing deck of cards. +would add a new folder into a file cabinet. It is a complex function, and we illustrate it below with a set of examples that sequentially reveal its features. @@ -326,7 +326,7 @@ \begin{itemize} \item Each row of the resulting data set represents a time interval (time1, time2] which is open on the left and closed on the right. - Covariate values for that row are the covariate values that apply + Covariate values for each row are the covariate values that apply over that interval. \item The event variable for each row $i$ is 1 if the time interval ends with an event and 0 otherwise. @@ -417,8 +417,8 @@ \subsection{Stanford heart transplant} The \code{jasa} data set contains information from the Stanford heart -transplant study, in the form that it appeared in the -paper of Crowley and Hu \cite{Crowley77}. +transplant study, in the form that it appeared in the J American Statistical +Association paper of Crowley and Hu \cite{Crowley77}. The data set has one line per subject which contains the baseline covariates along with dates of enrollment, transplant, and death or last follow-up. @@ -439,7 +439,7 @@ treatment for days 1--10 and as transplanted starting on day 11. That is, except for patient 38 who died on the same day as their procedure. They should be treated as a transplant death; the problem is resolved by - moving this transplant back .5 day. + moving this forward in time by .5 day. \item The treatment coefficients in table 6.1 of the definitive analysis found in @@ -497,7 +497,9 @@ \code{tdc} or \code{cumtdc} call that has two arguments and 0 in all other cases. If the variable being created is already a part of data1, -then our updates make changes to that variable. +then our updates make changes to that variable; in the case of +an \code{event} call it results in an addition, for \code{tdc} +it results in wholesale replacement of any prior values. Be careful of this. This feature is what allowed for the \code{infection} indicator to be build up incrementally in the cgd example above, but quite surprising results can occur when you