[med-svn] [r-cran-cmprsk] 04/06: New upstream version 2.2-7
Andreas Tille
tille at debian.org
Thu Oct 19 08:56:35 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-cmprsk.
commit c7696e2ac81032acdb03e2bd892d2a1f729d25e2
Author: Andreas Tille <tille at debian.org>
Date: Thu Oct 19 10:54:36 2017 +0200
New upstream version 2.2-7
---
COPYING | 340 +++++++++++
DESCRIPTION | 19 +
MD5 | 21 +
NAMESPACE | 19 +
R/cmprsk.R | 528 +++++++++++++++++
debian/README.test | 8 -
debian/changelog | 12 -
debian/compat | 1 -
debian/control | 25 -
debian/copyright | 30 -
debian/docs | 3 -
debian/rules | 3 -
debian/source/format | 1 -
debian/tests/control | 3 -
debian/tests/run-unit-test | 39 --
debian/upstream/metadata | 10 -
debian/watch | 2 -
man/crr.Rd | 167 ++++++
man/cuminc.Rd | 92 +++
man/extract.cuminc.Rd | 27 +
man/plot.cuminc.Rd | 73 +++
man/plot.predict.crr.Rd | 45 ++
man/predict.crr.Rd | 44 ++
man/print.crr.Rd | 28 +
man/print.cuminc.Rd | 30 +
man/summary.crr.Rd | 64 +++
man/timepoints.Rd | 34 ++
src/cincsub.f | 104 ++++
src/crr.f | 577 +++++++++++++++++++
src/crstm.f | 244 ++++++++
src/tpoi.f | 34 ++
tests/Rplots.ps | 1340 ++++++++++++++++++++++++++++++++++++++++++++
tests/test.R | 53 ++
tests/test.Rout.save | 476 ++++++++++++++++
34 files changed, 4359 insertions(+), 137 deletions(-)
diff --git a/COPYING b/COPYING
new file mode 100644
index 0000000..d60c31a
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,340 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Library General
+Public License instead of this License.
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..389637b
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,19 @@
+Package: cmprsk
+Version: 2.2-7
+Date: 2014-jun-17
+Title: Subdistribution Analysis of Competing Risks
+Author: Bob Gray <gray at jimmy.harvard.edu>
+Maintainer: Bob Gray <gray at jimmy.harvard.edu>
+Depends: R (>= 2.15.0), survival
+Description: Estimation, testing and regression modeling of
+ subdistribution functions in competing risks, as described in Gray
+ (1988), A class of K-sample tests for comparing the cumulative
+ incidence of a competing risk, Ann. Stat. 16:1141-1154, and Fine JP and
+ Gray RJ (1999), A proportional hazards model for the subdistribution
+ of a competing risk, JASA, 94:496-509.
+License: GPL (>= 2)
+URL: http://www.r-project.org
+Packaged: 2014-06-17 18:45:38 UTC; gray
+NeedsCompilation: yes
+Repository: CRAN
+Date/Publication: 2014-06-17 23:16:43
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..13d668c
--- /dev/null
+++ b/MD5
@@ -0,0 +1,21 @@
+94d55d512a9ba36caa9b7df079bae19f *COPYING
+024326158b9ee181904ce0688bb35ad0 *DESCRIPTION
+25e59606a24d9d76dd829a09d9a90c54 *NAMESPACE
+3e0e744d860b0ba3098fbcc8238e053e *R/cmprsk.R
+ab1f8f638eba7c7a13db3225977589a3 *man/crr.Rd
+8e35bb36130364553fcd0d80d18283c1 *man/cuminc.Rd
+9aa901885ef9cb0653c5ad154cb08735 *man/extract.cuminc.Rd
+6fa0c161e6dd9b2b9db731b70b2fa3bf *man/plot.cuminc.Rd
+d055318772f89bf1d6ce0c14d5d63331 *man/plot.predict.crr.Rd
+e1377154ce50530a986b51d422d50424 *man/predict.crr.Rd
+9eee30941b7ba829928c199fb50e16f6 *man/print.crr.Rd
+446fe5bacfb526891869147d40563bf6 *man/print.cuminc.Rd
+6fc86389f06d65ca888dc40e8961b225 *man/summary.crr.Rd
+3d843553233387c43d273a8f721722bb *man/timepoints.Rd
+c53028d318cc93d99254bb7be508b04c *src/cincsub.f
+abcc0f1aca5c56602513cfb1451bd19a *src/crr.f
+d02825e3238b85f763cea36899825634 *src/crstm.f
+a5c0126cbcc253811db0e486915e6ca9 *src/tpoi.f
+ffcf7be80bea500c4f61cc8bd37fee03 *tests/Rplots.ps
+8c8e6e952b759046af69be4d750bdc63 *tests/test.R
+63a0169768e776a0776a040d826c2116 *tests/test.Rout.save
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..c9c729e
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,19 @@
+# Export all names
+# exportPattern(".")
+useDynLib(cmprsk)
+export("crr","summary.crr","print.summary.crr","predict.crr",
+ "plot.predict.crr","print.crr","cuminc","print.cuminc","timepoints",
+ "plot.cuminc","[.cuminc")
+#export("crr","cuminc","timepoints")
+# Import all packages listed as Imports or Depends
+import(
+ survival
+)
+S3method(summary,crr)
+S3method(print,crr)
+S3method(print,summary.crr)
+S3method(predict,crr)
+S3method(print,cuminc)
+S3method(plot,cuminc)
+S3method('[',cuminc)
+S3method(plot,predict.crr)
diff --git a/R/cmprsk.R b/R/cmprsk.R
new file mode 100644
index 0000000..092c263
--- /dev/null
+++ b/R/cmprsk.R
@@ -0,0 +1,528 @@
+# Copyright (C) 2000 Robert Gray
+# distributed under the terms of the GNU public license
+crr <-
+# function for regression modeling of subdistribution functions
+# arguments:
+# ftime = vector of failure/censoring times
+# fstatus = vector with a unique code for each failure type and a
+# separate code for censored observations
+# cov1 = (nobs x ncovs) matrix of fixed covariates
+# cov2 = matrix of covariates multiplied by functions of time;
+# if used, often these covariates would also appear in cov1,
+# to give a prop hazards effect plus a time interaction
+# tf = functions of time. A function that takes a vector of times as
+# an argument and returns a matrix whose jth column is the value of
+# the time function corresponding to the jth column of cov2 evaluated
+# at the input time vector. At time tk, the
+# model includes the term cov2[,j]*tfs(tk)[,j] as a covariate.
+# cengroup = vector with different values for each group with
+# a distinct censoring distribution (the censoring distribution
+# is estimated separately within these groups)
+# failcode = code of fstatus that denotes the failure type of interest
+# cencode = code of fstatus that denotes censored observations
+# subset = logical vector length(ftime) indicating which cases to include
+# na.action = function defining action to take for cases that have NA for
+# any of ftime, fstatus, cov1, cov2 cengroup, or subset.
+# gtol = iteration stops when a function of the gradient is < gtol.
+# maxiter = maximum # of iterations in Newton algorithm (0 computes
+# scores and var at init, but performs no iterations)
+# init = initial values of regression parameters
+function(ftime,fstatus,cov1,cov2,tf,cengroup,failcode=1,cencode=0,
+ subset,na.action=na.omit,gtol=1e-6,maxiter=10,init,variance=TRUE) {
+ ## LS
+ call <- match.call()
+ cov1.name <- deparse(substitute(cov1))
+ cov1.vars <- cov2.vars <- NULL
+ if(!missing(cov1))
+ { cov1.vars <- colnames(as.matrix(cov1)) }
+ cov2.name <- deparse(substitute(cov2))
+ if(!missing(cov2))
+ { cov2.vars <- colnames(as.matrix(cov2)) }
+ ##
+ d <- data.frame(ftime=ftime,fstatus=fstatus,
+ cengroup=if (missing(cengroup)) rep(1,length(fstatus)) else cengroup)
+ if (!missing(cov1)) {
+ cov1 <- as.matrix(cov1)
+ nc1 <- ncol(cov1)
+ d <- cbind(d,cov1)
+ } else {nc1 <- 0}
+ if (!missing(cov2)) {
+ cov2 <- as.matrix(cov2)
+ nc2 <- ncol(cov2)
+ d <- cbind(d,cov2)
+ } else {nc2 <- 0}
+ if (!missing(subset)) d <- d[subset,]
+ tmp <- nrow(d)
+ d <- na.action(d)
+ nmis <- 0
+ if (nrow(d) != tmp) {
+ nmis <- tmp-nrow(d)
+ cat(format(nmis),'cases omitted due to missing values\n')
+ }
+ d <- d[order(d$ftime),]
+ ftime <- d$ftime
+ cenind <- ifelse(d$fstatus==cencode,1,0)
+ fstatus <- ifelse(d$fstatus==failcode,1,2*(1-cenind))
+ ucg <- sort(unique.default(d$cengroup))
+ cengroup <- match(d$cengroup,ucg)
+ ncg <- length(ucg)
+ uuu <- matrix(0,nrow=ncg,ncol=length(ftime))
+ for (k in 1:ncg) {
+ u <- do.call('survfit',list(formula=Surv(ftime,cenind)~1,data=
+ data.frame(ftime,cenind,cengroup),subset=cengroup==k))
+### note: want censring dist km at ftime-
+ u <- approx(c(0,u$time,max(u$time)*(1+10*.Machine$double.eps)),c(1,u$surv,
+ 0),xout=ftime*(1-100*.Machine$double.eps),method='constant',f=0,rule=2)
+ uuu[k,1:length(u$y)] <- u$y
+# u <- summary(u,times=sort(ftime*(1-.Machine$double.eps)))
+# uuu[k,1:length(u$surv)] <- u$surv
+ }
+ uft <- sort(unique(ftime[fstatus==1]))
+ ndf <- length(uft)
+ if (nc2 == 0) {
+ cov1 <- as.matrix(d[,(1:nc1)+3])
+ np <- nc1
+ npt <- 0
+ cov2 <- 0
+ tfs <- 0
+ } else if (nc1 == 0) {
+ cov2 <- as.matrix(d[,(1:nc2)+3+nc1])
+ npt <- np <- nc2
+ cov1 <- 0
+ tfs <- tf(uft)
+ } else {
+ cov1 <- as.matrix(d[,(1:nc1)+3])
+ cov2 <- as.matrix(d[,(1:nc2)+3+nc1])
+ npt <- nc2
+ np <- nc1+nc2
+ tfs <- tf(uft)
+ }
+### start of nr
+ if (missing(init)) b <- rep(0,np)
+ else b <- init
+ stepf <- .5
+ for (ll in 0:maxiter) {
+ z <- .Fortran('crrfsv',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(b),
+ double(1),double(np),double(np*np),double(np),double(np),
+ double(np*np),PACKAGE = "cmprsk")[15:17]
+ if (max(abs(z[[2]])*pmax(abs(b),1)) < max(abs(z[[1]]),1)*gtol) {
+ converge <- TRUE
+ break
+ }
+ if (ll==maxiter) {
+ converge <- FALSE
+ break
+ }
+ h <- z[[3]]
+ dim(h) <- c(np,np)
+### better to guarantee a pd factorization, but
+### matrix should be pd except in rare circumstances
+ sc <- -solve(h,z[[2]])
+ bn <- b+sc
+ fbn <- .Fortran('crrf',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(bn),
+ double(1),double(np),PACKAGE = "cmprsk")[[15]]
+# backtracking loop
+ i <- 0
+ while (is.na(fbn) || fbn>z[[1]]+(1e-4)*sum(sc*z[[2]])) {
+ i <- i+1
+ sc <- sc*stepf
+ bn <- b+sc
+ fbn <- .Fortran('crrf',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(bn),
+ double(1),double(np),PACKAGE = "cmprsk")[[15]]
+ if (i>20) break
+ }
+ if (i>20) {
+ converge <- FALSE
+ break
+ }
+ b <- c(bn)
+ }
+ if (variance) {
+ v <- .Fortran('crrvv',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(b),
+ double(np*np),double(np*np),double(np*np),
+ double(length(ftime)*(np+1)),double(np),double(np*ncg),
+ double(2*np),double(ncg*np),
+ integer(ncg),double(ncg*np),double(ncg),
+ PACKAGE = "cmprsk")[15:16]
+ dim(v[[2]]) <- dim(v[[1]]) <- c(np,np)
+ h0 <- v[[1]]
+ h <- solve(v[[1]])
+ v <- h %*% v[[2]] %*% t(h)
+ r <- .Fortran('crrsr',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(b),
+ double(ndf*np),double(np),double(np),PACKAGE = "cmprsk")[[15]]
+ r <- t(matrix(r,nrow=np))
+##
+ } else {
+ v <- h <- h0 <- matrix(NA,np,np)
+ r <- NULL
+# bj <- NULL
+ }
+ nobs <- length(ftime)
+ b0 <- rep(0,length(b))
+ fb0 <- .Fortran('crrf',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(b0),
+ double(1),double(np),PACKAGE = "cmprsk")[[15]]
+ bj <- .Fortran('crrfit',as.double(ftime),as.integer(fstatus),
+ as.integer(length(ftime)),as.double(cov1),as.integer(np-npt),
+ as.integer(np),as.double(cov2),as.integer(npt),
+ as.double(tfs),as.integer(ndf),as.double(uuu),
+ as.integer(ncg),as.integer(cengroup),as.double(b),
+ double(ndf),double(np),PACKAGE = "cmprsk")[[15]]
+ if (nc1>0) {
+ x1 <- paste(cov1.name, 1:nc1, sep="")
+ if(is.null(cov1.vars)) cov1.vars <- x1
+ else cov1.vars <- ifelse(cov1.vars=="",x1,cov1.vars)
+ }
+ if(nc2 > 0) {
+ x1 <- paste(cov2.name, 1:nc2, sep="")
+ if (is.null(cov2.vars)) cov2.vars <- x1
+ else cov2.vars <- ifelse(cov2.vars=="",x1,cov2.vars)
+ x1 <- paste('tf',1:nc2,sep='')
+ x2 <- colnames(tfs)
+ if (!is.null(x2)) x1 <- ifelse(x2=="",x1,x2)
+ cov2.vars <- paste(cov2.vars, x1, sep="*")
+ }
+ names(b) <- c(cov1.vars, cov2.vars)
+ ##
+ z <- list(coef=b,loglik=-z[[1]],score=-z[[2]],inf=h0,
+ var=v,res=r,uftime=uft,bfitj=bj,
+ tfs=as.matrix(tfs),converged=converge,call = call, n = nobs,
+ n.missing = nmis, loglik.null = -fb0,invinf=h)
+ class(z) <- 'crr'
+ z
+}
+
+"summary.crr" <- function(object, conf.int = 0.95, digits = max(options()$digits - 5, 2), ...)
+{
+ beta <- object$coef
+ se <- sqrt(diag(object$var))
+ out <- list(call = object$call, converged = object$converged,
+ n = object$n, n.missing = object$n.missing,
+ loglik = object$loglik)
+ tmp <- cbind(beta, exp(beta), se, beta/se,
+ signif(2 * (1 - pnorm(abs(beta)/se)), digits))
+ dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)",
+ "se(coef)", "z", "p-value"))
+ out$coef <- tmp
+ if(conf.int)
+ { a <- (1 - conf.int)/2
+ a <- c(a, 1 - a)
+ z <- qnorm(a)
+ tmp <- cbind(exp(beta), exp(-beta),
+ exp(beta + z[1] * se), exp(beta + z[2] * se))
+ dimnames(tmp) <- list(names(beta), c("exp(coef)", "exp(-coef)",
+ paste(format(100*a, trim = TRUE,
+ scientific = FALSE,
+ digits = 3), "%", sep="")))
+ out$conf.int <- tmp
+ }
+ df <- length(beta)
+ logtest <- -2 * (object$loglik.null - object$loglik)
+ out$logtest <- c(test = logtest, df = df)
+ # out$rsq <- c(rsq = 1 - exp(-logtest/object$n),
+ # maxrsq = 1 - exp(2 * object$loglik.null/object$n))
+ class(out) <- "summary.crr"
+ out
+}
+
+"print.summary.crr" <- function (x, digits = max(options()$digits - 4, 3), ...)
+{
+ cat("Competing Risks Regression\n\n")
+ if(!is.null(x$call))
+ { cat("Call:\n")
+ dput(x$call)
+ cat("\n")
+ }
+ if(!x$converged)
+ { cat("crr converged:", x$converged, "\n")
+ return()
+ }
+ savedig <- options(digits = digits)
+ on.exit(options(savedig))
+ print(x$coef)
+ cat("\n")
+ print(x$conf.int)
+ cat("\n")
+ cat("Num. cases =", x$n)
+ if(x$n.missing > 0)
+ cat(" (", x$n.missing, " cases omitted due to missing values)", sep="")
+ cat("\n")
+ # cat("Rsquare =", format(round(x$rsq["rsq"], 3)), " (max possible =",
+ # format(round(x$rsq["maxrsq"], 3)), ")\n")
+ cat("Pseudo Log-likelihood =", x$loglik, "\n")
+ cat("Pseudo likelihood ratio test = ", format(round(x$logtest["test"], 2)),
+ " on ", x$logtest["df"], " df,", "\n", sep = "")
+# " p-value = ", format(x$logtest["pvalue"]), "\n", sep = "")
+ invisible()
+}
+
+predict.crr <-
+# for a crr object x, estimates subdistributions at covariate
+# combinations given by rows of cov1 and cov2. The terms in cov1
+# cov2 must correspond exactly to the corresponding call to crr.
+ function(object,cov1,cov2,...) {
+ if (is.null(object$bfitj)) stop('predict requires variance=TRUE in crr')
+ np <- length(object$coef)
+ if (length(object$tfs)<=1) {
+ if (length(object$coef)==length(cov1)) lhat <- cumsum(exp(sum(cov1*object$coef))*object$bfitj)
+ else {
+ cov1 <- as.matrix(cov1)
+ lhat <- matrix(0,nrow=length(object$uftime),ncol=nrow(cov1))
+ for (j in 1:nrow(cov1)) lhat[,j] <- cumsum(exp(sum(cov1[j,]*object$coef))*object$bfitj)
+ }
+ } else {
+ if (length(object$coef)==ncol(as.matrix(object$tfs))) {
+ if (length(object$coef)==length(cov2))
+ lhat <- cumsum(exp(object$tfs %*% c(cov2*object$coef))*object$bfitj)
+ else {
+ cov2 <- as.matrix(cov2)
+ lhat <- matrix(0,nrow=length(object$uftime),ncol=nrow(cov1))
+ for (j in 1:nrow(cov2)) lhat[,j] <-
+ cumsum(exp(object$tfs %*% c(cov2[j,]*object$coef))*object$bfitj)
+ }
+ } else {
+ if (length(object$coef)==length(cov1)+length(cov2))
+ lhat <- cumsum(exp(sum(cov1*object$coef[1:length(cov1)])+object$tfs %*%
+ c(cov2*object$coef[(np-length(cov2)+1):np]))*object$bfitj)
+ else {
+ cov1 <- as.matrix(cov1)
+ cov2 <- as.matrix(cov2)
+ lhat <- matrix(0,nrow=length(object$uftime),ncol=nrow(cov1))
+ for (j in 1:nrow(cov1)) lhat[,j] <-
+ cumsum(exp(sum(cov1[j,]*object$coef[1:ncol(cov1)])+object$tfs %*%
+ c(cov2[j,]*object$coef[(np-ncol(cov2)+1):np]))*object$bfitj)
+ }
+ }
+ }
+ lhat <- cbind(object$uftime,1-exp(-lhat))
+ class(lhat) <- 'predict.crr'
+ lhat
+ }
+plot.predict.crr <-
+# plots estimated subdistributions from predict.crr
+ function(x,lty=1:(ncol(x)-1),color=1,ylim=c(0,max(x[,-1])),xmin=0,xmax=max(x[,1]),...) {
+ if (length(lty)<ncol(x)-1) lty <- rep(lty[1],ncol(x)-1)
+ if (length(color)<ncol(x)-1) color <- rep(color[1],ncol(x)-1)
+ if (xmax<max(x[,1])) x <- x[x[,1]<xmax,]
+ times <- c(xmin,rep(x[,1],rep(2,nrow(x))),xmax)
+ plot(c(xmin,xmax),ylim,type='n',...)
+ for (j in 2:ncol(x)) lines(times,c(0,0,rep(x[,j],rep(2,nrow(x)))),lty=lty[j-1],col=color[j-1])
+}
+print.crr <-
+# prints a summary of the crr fit x
+ function(x,...) {
+ cat('convergence: ',x$converged,'\n')
+ cat('coefficients:\n')
+ print(signif(x$coef,4),...)
+ v <- sqrt(diag(x$var))
+ cat('standard errors:\n')
+ print(signif(v,4),...)
+ v <- 2*(1-pnorm(abs(x$coef)/v))
+ cat('two-sided p-values:\n')
+ print(signif(v,2),...)
+ invisible()
+}
+
+cuminc <- function(ftime,fstatus,group,strata,rho=0,cencode=0,subset,na.action=na.omit) {
+# ftime=failure times, fstatus=variable which indicates the type
+# of failure (and cens), group is the group variable, strata=
+# strata variables for the tests (omit if none), rho is the
+# power of the weight function used in the tests, cencode is the value
+# of fstatus which indicates that a time is censored (default is 0)
+# subset = logical vector length(ftime) indicating which cases to include
+# na.action = function defining action to take for cases that have NA for
+# any of ftime, fstatus, group, strata, or subset.
+# output is a list giving the estmated cuminc
+# functions (times, function values, variances) for each group, and
+# a component
+# values of the test statistics for comparing each cause among the
+# groups. the tests are stratified (if strata specified).
+# check lengths, and status of group and strata
+ d <- data.frame(time=ftime,cause=fstatus,
+ group=as.factor(if (missing(group)) rep(1,length(ftime)) else group),
+ strata=as.factor(if (missing(strata)) rep(1,length(ftime)) else strata))
+ if (!missing(subset)) d <- d[subset,]
+ tmp <- nrow(d)
+ d <- na.action(d)
+ if (nrow(d) != tmp) cat(format(tmp-nrow(d)),'cases omitted due to missing values\n')
+ no <- nrow(d)
+ cg <- " "
+ nst <- length(levels(d$strata))
+ d <- d[order(d$time),]
+ ugg <- table(d$group)
+ d$group <- factor(d$group,names(ugg)[ugg>0])
+ ugg <- levels(d$group)
+ censind <- ifelse(d$cause==cencode,0,1)
+ uc <- table(d$cause[censind==1])
+ if (is.factor(d$cause)) uclab <- names(uc)[uc>0]
+ else uclab <- as.numeric(names(uc)[uc>0])
+ nc <- length(uclab)
+ ng <- length(ugg)
+ if (ng>1) {
+ ng1 <- ng-1
+ ng2 <- ng*ng1/2
+ v <- matrix(0,nrow=ng1,ncol=ng1)
+ storage.mode(v) <- "double"
+ vt <- double(ng2)
+ s <- double(ng1)
+ }
+ pf <- vector("list",ng*nc)
+ stat <- double(nc)
+ l <- 0
+ for (ii in 1:nc) {
+ causeind <- ifelse(d$cause==uclab[ii],1,0)
+ for (jj in 1:length(ugg)) {
+ cg <- c(cg,paste(ugg[jj],uclab[ii]))
+ l <- l+1
+ cgind <- d$group==ugg[jj]
+ ncg <- length(cgind[cgind])
+ n2 <- length(unique(d$time[cgind & causeind==1]))
+ n2 <- 2*n2+2
+ tmp <- double(n2)
+ z <- .Fortran("cinc",as.double(d$time[cgind]),as.integer(censind[cgind]),
+ as.integer(causeind[cgind]),as.integer(ncg),
+ x=tmp,f=tmp,v=tmp,PACKAGE = "cmprsk")
+ pf[[l]] <- list(time=z$x,est=z$f,var=z$v)
+ }
+ if (ng>1) {
+ causeind <- 2*censind-causeind
+ z2 <- .Fortran("crstm",as.double(d$time),as.integer(causeind),
+ as.integer(d$group),as.integer(d$strata),as.integer(no),
+ as.double(rho),as.integer(nst),as.integer(ng),s,v,
+ as.double(d$time),as.integer(causeind),as.integer(d$group),
+ vt,s,vt,double((4+3*ng)*ng),integer(4*ng),PACKAGE = "cmprsk")
+ stat[ii] <- -1
+ a <- qr(z2[[10]])
+ if (a$rank==ncol(a$qr)) {
+ b <- diag(dim(a$qr)[1])
+ stat[ii] <- z2[[9]]%*%qr.coef(a,b)%*%z2[[9]]
+ }
+ }
+ }
+ names(pf) <- cg[2:length(cg)]
+ if (ng>1) {
+ names(stat) <- uclab
+ stat <- list(Tests=cbind(stat=stat,pv=1-pchisq(stat,ng-1),df=rep(ng-1,length(stat))))
+ pf <- c(pf,stat)
+ }
+ attr(pf, "class") <- "cuminc"
+ pf
+}
+
+print.cuminc <- function(x,ntp=4,maxtime,...) {
+ if (!is.null(x$Tests)) {
+ cat('Tests:\n')
+ print(x$Tests)
+ nc <- length(x)-1
+ } else {
+ nc <- length(x)
+ }
+ if (missing(maxtime)) {
+ maxtime <- 0
+ for (i in 1:nc) maxtime <- max(maxtime,x[[i]]$time)
+ }
+ tp <- pretty(c(0,maxtime),ntp+1)
+ cat('Estimates and Variances:\n')
+ print(timepoints(x,tp[-c(1,length(tp))]),...)
+ invisible()
+}
+
+timepoints <- function(w,times) {
+# w=list (see cuminc or km), times= times you want estimates for.
+# output is a list with components est giving the estimates, var giving
+# the variances,
+ if (!is.null(w$Tests)) w <- w[names(w) != 'Tests']
+ ng <- length(w)
+ times <- sort(unique(times))
+ nt <- length(times)
+ storage.mode(times) <- "double"
+ storage.mode(nt) <- "integer"
+ ind <- matrix(0,ncol=nt,nrow=ng)
+ oute <- matrix(NA,ncol=nt,nrow=ng)
+ outv <- oute
+ storage.mode(ind) <- "integer"
+ slct <- rep(TRUE,ng)
+ for (i in 1:ng) {
+ if (is.null((w[[i]])$est)) { slct[i] <- FALSE} else {
+ z <- .Fortran("tpoi",as.double(w[[i]][[1]]),
+ as.integer(length(w[[i]][[1]])),ind[i,],times,nt,PACKAGE = "cmprsk")
+ ind[i,] <- z[[3]]
+ oute[i,ind[i,]>0] <- w[[i]][[2]][z[[3]]]
+ if (length(w[[i]])>2) outv[i,ind[i,]>0] <- w[[i]][[3]][z[[3]]]
+ }
+ }
+ dimnames(oute) <- list(names(w)[1:ng],as.character(times))
+ dimnames(outv) <- dimnames(oute)
+ list(est=oute[slct,,drop=FALSE],var=outv[slct,,drop=FALSE])
+}
+
+plot.cuminc <- function(x,main=" ",curvlab,ylim=c(0,1),xlim,wh=2,xlab="Years",
+ylab="Probability",lty=1:length(x),color=1,lwd = par('lwd'),...) {
+# x is a list containing curves to be plotted. Each component of
+# x is a list with the first component containing the x values
+# and the second component the y values. main = main title in the plot
+# curvlab=curve labels (vector), wh=where curve labels are plotted
+# 1=lower left 2=upper left 3=upper right 4=lower right
+ if (!is.null(x$Tests)) x <- x[names(x) != 'Tests']
+ nc <- length(x)
+ if (length(lty) < nc) lty <- rep(lty[1],nc) else lty <- lty[1:nc]
+ if (length(lwd) < nc) lwd <- rep(lwd[1],nc) else lwd <- lwd[1:nc]
+ if (length(color) < nc) color <- rep(color[1],nc) else color <- color[1:nc]
+ if (missing(curvlab)) {
+ if (mode(names(x))=="NULL") {
+ curvlab <- as.character(1:nc) }
+ else curvlab <- names(x)[1:nc]
+ }
+ if (missing(xlim)) {
+ xmax <- 0
+ for (i in 1:nc) {
+ xmax <- max(c(xmax,x[[i]][[1]]))
+ }
+ xlim <- c(0,xmax)
+ }
+ plot(x[[1]][[1]],x[[1]][[2]],type="n",ylim=ylim,xlim=xlim,
+ main=main,xlab=xlab,ylab=ylab,bty="l",...)
+ if (length(wh) != 2) {
+ wh <- c(xlim[1],ylim[2])
+ }
+ u <- list(...)
+ if (length(u)>0) {
+ i <- pmatch(names(u),names(formals(legend)),0)
+ do.call('legend',c(list(x=wh[1],y=wh[2],legend=curvlab,col=color,lty=lty,lwd=lwd,bty="n",bg=-999999),u[i>0]))
+ } else {
+ do.call('legend',list(x=wh[1],y=wh[2],legend=curvlab,col=color,lty=lty,lwd=lwd,bty="n",bg=-999999))
+ }
+# legend(wh[1],wh[2],legend=curvlab,col=color,lty=lty,bty="n",bg=-999999,...)
+ for (i in 1:nc) {
+ lines(x[[i]][[1]],x[[i]][[2]],lty=lty[i],col=color[i],lwd=lwd[i],...)
+ }
+}
+
+"[.cuminc" <- function(x,i,...) {
+ x <- NextMethod("[")
+ class(x) <- 'cuminc'
+ x
+}
diff --git a/debian/README.test b/debian/README.test
deleted file mode 100644
index 55a9142..0000000
--- a/debian/README.test
+++ /dev/null
@@ -1,8 +0,0 @@
-Notes on how this package can be tested.
-────────────────────────────────────────
-
-To run the unit tests provided by the package you can do
-
- sh run-unit-test
-
-in this directory.
diff --git a/debian/changelog b/debian/changelog
deleted file mode 100644
index 1ded39f..0000000
--- a/debian/changelog
+++ /dev/null
@@ -1,12 +0,0 @@
-r-cran-cmprsk (2.2-7-2) unstable; urgency=medium
-
- * Make test more tolerant against different output to compare with
- * cme fix dpkg-control
-
- -- Andreas Tille <tille at debian.org> Fri, 29 Apr 2016 09:10:40 +0200
-
-r-cran-cmprsk (2.2-7-1) unstable; urgency=low
-
- * Initial release (closes: #799920)
-
- -- Andreas Tille <tille at debian.org> Thu, 24 Sep 2015 11:22:20 +0200
diff --git a/debian/compat b/debian/compat
deleted file mode 100644
index ec63514..0000000
--- a/debian/compat
+++ /dev/null
@@ -1 +0,0 @@
-9
diff --git a/debian/control b/debian/control
deleted file mode 100644
index 6f4f7ac..0000000
--- a/debian/control
+++ /dev/null
@@ -1,25 +0,0 @@
-Source: r-cran-cmprsk
-Maintainer: Debian Med Packaging Team <debian-med-packaging at lists.alioth.debian.org>
-Uploaders: Andreas Tille <tille at debian.org>
-Section: gnu-r
-Testsuite: autopkgtest
-Priority: optional
-Build-Depends: debhelper (>= 9),
- cdbs,
- r-base-dev,
- r-cran-survival
-Standards-Version: 3.9.8
-Vcs-Browser: https://anonscm.debian.org/viewvc/debian-med/trunk/packages/R/r-cran-cmprsk/trunk/
-Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/R/r-cran-cmprsk/trunk/
-Homepage: https://cran.r-project.org/web/packages/cmprsk/
-
-Package: r-cran-cmprsk
-Architecture: any
-Depends: ${shlibs:Depends},
- ${R:Depends},
- r-cran-survival
-Description: GNU R subdistribution analysis of competing risks
- This GNU R package supports estimation, testing and regression modeling
- of subdistribution functions in competing risks, as described in Gray
- (1988), A class of K-sample tests for comparing the cumulative incidence
- of a competing risk.
diff --git a/debian/copyright b/debian/copyright
deleted file mode 100644
index b78c98f..0000000
--- a/debian/copyright
+++ /dev/null
@@ -1,30 +0,0 @@
-Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: cmprsk
-Upstream-Contact: Bob Gray <gray at jimmy.harvard.edu>
-Source: http://cran.r-project.org/src/contrib/
-
-Files: *
-Copyright: 1998-2014 Bob Gray <gray at jimmy.harvard.edu>
-License: GPL-2+
-
-Files: debian/*
-Copyright: 2015 Andreas Tille <tille at debian.org>
-License: GPL-2+
-
-License: GPL-2+
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 2 of the License, or
- (at your option) any later version.
- .
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- .
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
- .
- On Debian systems, the complete text of the GNU General Public
- License can be found in `/usr/share/common-licenses/GPL'.
-
diff --git a/debian/docs b/debian/docs
deleted file mode 100644
index 3adf0d6..0000000
--- a/debian/docs
+++ /dev/null
@@ -1,3 +0,0 @@
-debian/README.test
-debian/tests/run-unit-test
-tests
diff --git a/debian/rules b/debian/rules
deleted file mode 100755
index 2fbba2d..0000000
--- a/debian/rules
+++ /dev/null
@@ -1,3 +0,0 @@
-#!/usr/bin/make -f
-
-include /usr/share/R/debian/r-cran.mk
diff --git a/debian/source/format b/debian/source/format
deleted file mode 100644
index 163aaf8..0000000
--- a/debian/source/format
+++ /dev/null
@@ -1 +0,0 @@
-3.0 (quilt)
diff --git a/debian/tests/control b/debian/tests/control
deleted file mode 100644
index d2aa55a..0000000
--- a/debian/tests/control
+++ /dev/null
@@ -1,3 +0,0 @@
-Tests: run-unit-test
-Depends: @
-Restrictions: allow-stderr
diff --git a/debian/tests/run-unit-test b/debian/tests/run-unit-test
deleted file mode 100644
index 24dd25b..0000000
--- a/debian/tests/run-unit-test
+++ /dev/null
@@ -1,39 +0,0 @@
-#!/bin/sh -e
-
-pkg=r-cran-cmprsk
-
-# The saved result files do contain some differences in metadata and we also
-# need to ignore version differences of R
-filter() {
- grep -v -e '^R version' \
- -e '^Copyright (C)' \
- -e '^Platform: ' \
- -e '^ISBN 3' \
- -e '[Bb]uggy version of Kinderman-Ramage generator use' \
- -e '^Loading required package: splines' \
- $1 | \
- sed '/^> *proc\.time()$/,$d'
-}
-
-
-if [ "$ADTTMP" = "" ] ; then
- ADTTMP=`mktemp -d /tmp/${pkg}-test.XXXXXX`
-fi
-cd $ADTTMP
-cp /usr/share/doc/${pkg}/tests/* $ADTTMP
-find . -name "*.gz" -exec gunzip \{\} \;
-for htest in `ls *.R | sed 's/\.R$//'` ; do
- LC_ALL=C R --no-save < ${htest}.R 2>&1 | tee > ${htest}.Rout
- filter ${htest}.Rout.save > ${htest}.Rout.save_
- filter ${htest}.Rout > ${htest}.Rout_
- diff -u ${htest}.Rout.save_ ${htest}.Rout_
- if [ ! $? ] ; then
- echo "Test ${htest} failed"
- exit 1
- else
- echo "Test ${htest} passed"
- fi
-done
-rm -f $ADTTMP/*
-
-exit 0
diff --git a/debian/upstream/metadata b/debian/upstream/metadata
deleted file mode 100644
index 1c55833..0000000
--- a/debian/upstream/metadata
+++ /dev/null
@@ -1,10 +0,0 @@
-Reference:
- Author: Jason P. Fine and Robert J. Gray
- Title: "A proportional hazards model for the subdistribution of a competing risk"
- Journal: J Am Stat Assoc
- Year: 1999
- Volume: 94
- Number: 446
- Pages: 496-509
- DOI: 10.2307/2670170
- URL: http://www.tandfonline.com/doi/abs/10.1080/01621459.1999.10474144
diff --git a/debian/watch b/debian/watch
deleted file mode 100644
index 0e59700..0000000
--- a/debian/watch
+++ /dev/null
@@ -1,2 +0,0 @@
-version=3
-http://cran.r-project.org/src/contrib/cmprsk_([-\d.]*)\.tar\.gz
diff --git a/man/crr.Rd b/man/crr.Rd
new file mode 100644
index 0000000..b920fb4
--- /dev/null
+++ b/man/crr.Rd
@@ -0,0 +1,167 @@
+\name{crr}
+\alias{crr}
+\title{
+Competing Risks Regression
+}
+\description{
+regression modeling of subdistribution functions in competing risks
+}
+\usage{
+crr(ftime, fstatus, cov1, cov2, tf, cengroup, failcode=1, cencode=0,
+subset, na.action=na.omit, gtol=1e-06, maxiter=10, init, variance=TRUE)
+}
+\arguments{
+\item{ftime}{
+vector of failure/censoring times
+}
+\item{fstatus}{
+vector with a unique code for each failure type and a separate code for
+censored observations
+}
+\item{cov1}{
+matrix (nobs x ncovs) of fixed covariates (either cov1, cov2, or both
+are required)
+}
+\item{cov2}{
+matrix of covariates that will be multiplied by functions of time;
+if used, often these covariates would also appear in cov1
+to give a prop hazards effect plus a time interaction
+}
+\item{tf}{
+functions of time. A function that takes a vector of times as
+an argument and returns a matrix whose jth column is the value of
+the time function corresponding to the jth column of cov2 evaluated
+at the input time vector. At time \code{tk}, the
+model includes the term \code{cov2[,j]*tf(tk)[,j]} as a covariate.
+}
+\item{cengroup}{
+vector with different values for each group with
+a distinct censoring distribution (the censoring distribution
+is estimated separately within these groups). All data in one group, if
+missing.
+}
+\item{failcode}{
+code of fstatus that denotes the failure type of interest
+}
+\item{cencode}{
+code of fstatus that denotes censored observations
+}
+\item{subset}{
+ a logical vector specifying a subset of cases to include in the
+ analysis
+}
+\item{na.action}{
+ a function specifying the action to take for any cases missing any of
+ ftime, fstatus, cov1, cov2, cengroup, or subset.
+}
+\item{gtol}{
+ iteration stops when a function of the gradient is \code{< gtol}
+}
+\item{maxiter}{
+maximum number of iterations in Newton algorithm (0 computes
+scores and var at \code{init}, but performs no iterations)
+}
+\item{init}{
+ initial values of regression parameters (default=all 0)
+}
+\item{variance}{
+ If \code{FALSE}, then suppresses computation of the variance estimate
+ and residuals
+}
+}
+\value{
+ Returns a list of class crr, with components
+ \item{$coef}{the estimated regression coefficients}
+ \item{$loglik}{log pseudo-liklihood evaluated at \code{coef}}
+\item{$score}{derivitives of the log pseudo-likelihood evaluated at \code{coef}}
+\item{$inf}{-second derivatives of the log pseudo-likelihood}
+\item{$var}{estimated variance covariance matrix of coef}
+\item{$res}{matrix of residuals giving the
+contribution to each score (columns) at each unique failure time
+(rows)}
+\item{$uftime}{vector of unique failure times}
+\item{$bfitj}{jumps in the Breslow-type estimate of the underlying
+ sub-distribution cumulative hazard (used by predict.crr())}
+\item{$tfs}{the tfs matrix (output of tf(), if used)}
+\item{$converged}{TRUE if the iterative algorithm converged}
+\item{$call}{The call to crr}
+\item{$n}{The number of observations used in fitting the model}
+\item{$n.missing}{The number of observations removed from the input data
+ due to missing values}
+\item{$loglik.null}{The value of the log pseudo-likelihood when all the
+ coefficients are 0}
+\item{$invinf}{- inverse of second derivative matrix of the log pseudo-likelihood}
+}
+\details{
+Fits the 'proportional subdistribution hazards' regression model
+described in Fine and Gray (1999). This model directly assesses the
+effect of covariates on the subdistribution of a particular type of
+failure in a competing risks setting. The method implemented here is
+described in the paper as the weighted estimating equation.
+
+While the use of model formulas is not supported, the
+\code{model.matrix} function can be used to generate suitable matrices
+of covariates from factors, eg
+\code{model.matrix(~factor1+factor2)[,-1]} will generate the variables
+for the factor coding of the factors \code{factor1} and \code{factor2}.
+The final \code{[,-1]} removes the constant term from the output of
+\code{model.matrix}.
+
+The basic model assumes the subdistribution with covariates z is a
+constant shift on the complementary log log scale from a baseline
+subdistribution function. This can be generalized by including
+interactions of z with functions of time to allow the magnitude of the
+shift to change with follow-up time, through the cov2 and tfs
+arguments. For example, if z is a vector of covariate values, and uft
+is a vector containing the unique failure times for failures of the
+type of interest (sorted in ascending order), then the coefficients a,
+b and c in the quadratic (in time) model
+\eqn{az+bzt+zt^2}{a*z+b*z*t+c*z*t*t} can be fit
+by specifying \code{cov1=z}, \code{cov2=cbind(z,z)},
+\code{tf=function(uft) cbind(uft,uft*uft)}.
+
+This function uses an estimate of the survivor function of the
+censoring distribution to reweight contributions to the risk sets for
+failures from competing causes. In a generalization of the methodology
+in the paper, the censoring distribution can be estimated separately
+within strata defined by the cengroup argument. If the censoring
+distribution is different within groups defined by covariates in the
+model, then validity of the method requires using separate estimates of
+the censoring distribution within those groups.
+
+The residuals returned are analogous to the Schoenfeld residuals in
+ordinary survival models. Plotting the jth column of res against the
+vector of unique failure times checks for lack of fit over time in the
+corresponding covariate (column of cov1).
+
+If \code{variance=FALSE}, then %\code{predict.crr} cannot be used and
+some of the functionality in \code{summary.crr} and \code{print.crr}
+will be lost. This option can be useful in situations where crr is
+called repeatedly for point estimates, but standard errors are not
+required, such as in some approaches to stepwise model selection.
+}
+\references{
+Fine JP and Gray RJ (1999) A proportional hazards model for the
+subdistribution of a competing risk. JASA 94:496-509.
+}
+\seealso{
+\code{\link{predict.crr}} \code{\link{print.crr}} \code{\link{plot.predict.crr}}
+\code{\link{summary.crr}}
+}
+\examples{
+# simulated data to test
+set.seed(10)
+ftime <- rexp(200)
+fstatus <- sample(0:2,200,replace=TRUE)
+cov <- matrix(runif(600),nrow=200)
+dimnames(cov)[[2]] <- c('x1','x2','x3')
+print(z <- crr(ftime,fstatus,cov))
+summary(z)
+z.p <- predict(z,rbind(c(.1,.5,.8),c(.1,.5,.2)))
+plot(z.p,lty=1,color=2:3)
+crr(ftime,fstatus,cov,failcode=2)
+# quadratic in time for first cov
+crr(ftime,fstatus,cov,cbind(cov[,1],cov[,1]),function(Uft) cbind(Uft,Uft^2))
+#additional examples in test.R
+}
+\keyword{survival}
diff --git a/man/cuminc.Rd b/man/cuminc.Rd
new file mode 100644
index 0000000..bba5912
--- /dev/null
+++ b/man/cuminc.Rd
@@ -0,0 +1,92 @@
+\name{cuminc}
+\alias{cuminc}
+\title{
+Cumulative Incidence Analysis
+}
+\description{
+Estimate cumulative incidence functions from competing risks
+data and test equality across groups
+}
+\usage{
+cuminc(ftime, fstatus, group, strata, rho=0, cencode=0,
+subset, na.action=na.omit)
+}
+\arguments{
+\item{ftime}{
+failure time variable
+}
+\item{fstatus}{
+variable with distinct codes for different causes of failure
+and also a distinct code for censored observations
+}
+\item{group}{
+estimates will calculated within groups given by distinct values of this
+variable. Tests will compare these groups. If missing then treated as all
+one group (no test statistics)
+}
+\item{strata}{
+stratification variable. Has no effect on estimates. Tests will be
+stratified on this variable. (all data in 1 stratum, if missing)
+}
+\item{rho}{
+Power of the weight function used in the tests.
+}
+\item{cencode}{
+value of fstatus variable which indicates the failure time is censored.
+}
+\item{subset}{
+ a logical vector specifying a subset of cases to include in the
+ analysis
+}
+\item{na.action}{
+ a function specifying the action to take for any cases missing any of
+ ftime, fstatus, group, strata, or subset.
+}
+}
+\value{
+A list with components giving the subdistribution estimates for each
+cause in each group, and a component \code{Tests} giving the test
+statistics and p-values for comparing the subdistribution for each cause
+across groups (if the
+number of groups is \eqn{>}{>}1). The components giving the estimates
+have names that are a combination
+of the group name and the cause code.
+These components are also lists, with components
+\item{\code{time}}{ the times
+where the estimates are calculated}
+\item{\code{est}}{the estimated
+sub-distribution functions. These are step functions (all corners
+of the steps given), so they can be plotted using ordinary lines() commands.
+Estimates at particular times can be located using the timepoints()
+function.}
+\item{\code{var}}{the estimated variance of
+ the estimates, which are estimates of the asymptotic
+ variance of Aalen (1978). }
+}
+\references{
+Gray RJ (1988) A class of K-sample tests for comparing the cumulative
+incidence of a competing risk, ANNALS OF STATISTICS, 16:1141-1154.
+
+
+Kalbfleisch and Prentice (1980) THE ANALYSIS OF FAILURE TIME DATA, p 168-9.
+
+
+Aalen, O. (1978) Nonparametric estimation of partial transition
+probabilities in multiple decrement models, ANNALS OF STATISTICS,
+6:534-545.
+}
+\author{Robert Gray}
+\seealso{
+\code{\link{plot.cuminc}} \code{\link{timepoints}} \code{\link{print.cuminc}}
+}
+\examples{
+set.seed(2)
+ss <- rexp(100)
+gg <- factor(sample(1:3,100,replace=TRUE),1:3,c('a','b','c'))
+cc <- sample(0:2,100,replace=TRUE)
+strt <- sample(1:2,100,replace=TRUE)
+print(xx <- cuminc(ss,cc,gg,strt))
+plot(xx,lty=1,color=1:6)
+# see also test.R, test.out
+}
+\keyword{survival}
diff --git a/man/extract.cuminc.Rd b/man/extract.cuminc.Rd
new file mode 100644
index 0000000..1bcafd8
--- /dev/null
+++ b/man/extract.cuminc.Rd
@@ -0,0 +1,27 @@
+\name{[.cuminc}
+\alias{[.cuminc}
+\title{
+Subset method for lists of class cuminc
+}
+\description{
+A subset method that preserves the class of objects of class cuminc,
+allowing a subset of the curves to be selected.
+}
+\usage{
+\method{[}{cuminc}(x,i,\ldots)
+}
+\arguments{
+ \item{x}{object of class cuminc}
+ \item{i}{elements to extract}
+ \item{...}{not used}
+}
+
+\value{
+A list with selected components of \code{x}, with the class set to
+cuminc so cuminc methods can be applied.
+}
+\seealso{
+\code{\link{cuminc}} \code{\link{plot.cuminc}} \code{\link{print.cuminc}}
+}
+
+\keyword{survival}
diff --git a/man/plot.cuminc.Rd b/man/plot.cuminc.Rd
new file mode 100644
index 0000000..5c5344d
--- /dev/null
+++ b/man/plot.cuminc.Rd
@@ -0,0 +1,73 @@
+\name{plot.cuminc}
+\alias{plot.cuminc}
+\title{
+Create Labeled Cumulative Incidence Plots
+}
+\description{
+Plot method for cuminc. Creates labeled line plots from appropriate
+list input, for example, the output from \code{cuminc()}.
+}
+\usage{
+\method{plot}{cuminc}(x, main=" ", curvlab, ylim=c(0, 1), xlim, wh=2,
+xlab="Years", ylab="Probability", lty=1:length(x), color=1, lwd=par('lwd'),
+\dots)
+}
+\arguments{
+\item{x}{
+a list, with each component representing one curve in the plot. Each
+component of \code{x} is itself a list whose first component gives the x values
+and 2nd component the y values to be plotted. Although written for
+cumulative incidence curves, can in principle be used for any set of lines.
+}
+\item{main}{
+the main title for the plot.
+}
+\item{curvlab}{
+Curve labels for the plot. Default is \code{names(x)}, or if that is missing,
+\code{1:nc}, where \code{nc} is the number of curves in \code{x}.
+}
+\item{ylim}{
+yaxis limits for plot
+}
+\item{xlim}{
+xaxis limits for plot (default is 0 to the largest time in any of the
+curves)
+}
+\item{wh}{
+ if a vector of length 2, then the upper right coordinates of the
+ legend; otherwise the legend is placed in the upper right corner of
+ the plot
+}
+\item{xlab}{
+X axis label
+}
+\item{ylab}{
+ y axis label
+}
+\item{lty}{
+vector of line types. Default \code{1:nc} (\code{nc} is the number of
+curves in \code{x}). For color displays, \code{lty=1, color=1:nc}, might
+be more appropriate. If \code{length(lty)<nc}, then \code{lty[1]} is
+used for all.
+}
+\item{color}{
+ vector of colors. If \code{length(color)<nc}, then the \code{color[1]} is
+ used for all.
+ }
+\item{lwd}{
+ vector of line widths. If \code{length(lwd)<nc}, then \code{lwd[1]}
+ is used for all.
+ }
+ \item{...}{
+additional arguments passed to the initial call of the plot function.
+}}
+\value{
+No value is returned.
+}
+\seealso{ \code{\link{cuminc}} }
+%\examples{
+%#see help(cuminc)
+%}
+\keyword{survival}
+\keyword{hplot}
+% Converted by Sd2Rd version 1.10.
diff --git a/man/plot.predict.crr.Rd b/man/plot.predict.crr.Rd
new file mode 100644
index 0000000..754548f
--- /dev/null
+++ b/man/plot.predict.crr.Rd
@@ -0,0 +1,45 @@
+\name{plot.predict.crr}
+\alias{plot.predict.crr}
+\title{
+Plot estimated subdistribution functions
+}
+\description{
+plot method for \code{predict.crr}
+}
+\usage{
+\method{plot}{predict.crr}(x, lty=1:(ncol(x)-1), color=1,
+ylim=c(0, max(x[, -1])), xmin=0, xmax=max(x[, 1]), \dots)
+}
+\arguments{
+\item{x}{
+Output from \code{predict.crr}
+}
+\item{lty}{
+vector of line types. If length is \eqn{<}{<} \# curves, then
+\code{lty[1]} is used for all.
+}
+\item{color}{
+vector of line colors. If length is \eqn{<}{<} \# curves, then
+\code{color[1]} is used for all.
+}
+\item{ylim}{
+range of y-axis (vector of length two)
+}
+\item{xmin}{
+lower limit of x-axis (often 0, the default)
+}
+\item{xmax}{
+upper limit of x-axis
+}
+\item{...}{
+Other arguments to plot
+}}
+\section{Side Effects}{
+plots the subdistribution functions estimated by \code{predict.crr}, by
+default using a different line type for each curve
+}
+\seealso{
+\code{\link{crr}} \code{\link{predict.crr}}
+}
+\keyword{survival}
+% Converted by Sd2Rd version 1.10.
diff --git a/man/predict.crr.Rd b/man/predict.crr.Rd
new file mode 100644
index 0000000..2d3f9af
--- /dev/null
+++ b/man/predict.crr.Rd
@@ -0,0 +1,44 @@
+\name{predict.crr}
+\alias{predict.crr}
+\title{
+Estimate subdistribution functions from crr output
+}
+\description{
+predict method for crr
+}
+\usage{
+\method{predict}{crr}(object, cov1, cov2, \dots)
+}
+\arguments{
+\item{object}{
+output from crr
+}
+\item{cov1, cov2}{
+each row of cov1 and cov2 is a set of covariate values where the
+subdistribution should be estimated. The columns of cov1 and cov2 must
+be in the same order as in the original call to crr. Each must be
+given if present in the original call to crr.
+}
+\item{...}{
+additional arguments are ignored (included for compatibility with generic).
+}
+}
+\value{
+Returns a matrix with the unique type 1 failure times in the first
+column, and the other columns giving the estimated subdistribution
+function corresponding to the covariate combinations in the rows of cov1
+and cov2, at each failure time (the value that the estimate jumps to at
+that failure time).
+}
+\details{
+Computes \eqn{1-\exp(-B(t))}{1-exp(-B(t))}, where \eqn{B(t)}{B(t)} is
+the estimated cumulative
+sub-distribution hazard obtained for the specified covariate values,
+obtained from the Breslow-type estimate of the underlying hazard and
+the estimated regression coefficients.
+}
+\seealso{
+\code{\link{crr}} \code{\link{plot.predict.crr}}
+}
+\keyword{survival}
+% Converted by Sd2Rd version 1.10.
diff --git a/man/print.crr.Rd b/man/print.crr.Rd
new file mode 100644
index 0000000..3e7c34d
--- /dev/null
+++ b/man/print.crr.Rd
@@ -0,0 +1,28 @@
+\name{print.crr}
+\alias{print.crr}
+\title{
+prints summary of a crr object
+}
+\description{
+print method for crr objects
+}
+\usage{
+\method{print}{crr}(x, \dots)
+}
+\arguments{
+\item{x}{
+crr object (output from \code{crr()})
+}
+\item{...}{additional arguments to \code{print()}}
+}
+\details{
+prints the convergence status, the estimated coefficients, the
+estimated standard errors, and the two-sided p-values for the test of
+the individual coefficients equal to 0. (If convergence is false
+everything else may be meaningless.)
+}
+\seealso{
+\code{\link{crr}}
+}
+\keyword{survival}
+% Converted by Sd2Rd version 1.10.
diff --git a/man/print.cuminc.Rd b/man/print.cuminc.Rd
new file mode 100644
index 0000000..dbe2bd0
--- /dev/null
+++ b/man/print.cuminc.Rd
@@ -0,0 +1,30 @@
+\name{print.cuminc}
+\alias{print.cuminc}
+\title{Print cuminc objects}
+\description{
+A print method for objects of class cuminc (output from \code{cuminc()}).
+}
+\usage{
+\method{print}{cuminc}(x, ntp=4, maxtime, \dots)
+}
+\arguments{
+ \item{x}{an object of class cuminc}
+ \item{ntp}{number of timepoints where estimates are printed}
+ \item{maxtime}{the maximum timepoint where values are printed. The
+ default is the maximum time in the curves in \code{x}}
+ \item{...}{additional arguments to \code{print()}}
+}
+\details{
+Prints the test statistics and p-values (if present in \code{x}), and for each
+estimated cumulative incidence curve prints its value and estimated
+variance at a vector of times. The times are chosen between 0 and
+maxtime using the \code{pretty()} function.
+}
+\author{Robert Gray}
+
+\seealso{ \code{\link{cuminc}} }
+
+%\examples{
+%#see help(cuminc)
+%}
+\keyword{survival }%-- one or more ...
diff --git a/man/summary.crr.Rd b/man/summary.crr.Rd
new file mode 100644
index 0000000..8075273
--- /dev/null
+++ b/man/summary.crr.Rd
@@ -0,0 +1,64 @@
+\name{summary.crr}
+\alias{summary.crr}
+\alias{print.summary.crr}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Summary method for crr}
+\description{
+Generate and print summaries of crr output
+}
+\usage{
+\method{summary}{crr}(object, conf.int = 0.95, digits =
+max(options()$digits - 5, 2), ...)
+
+\method{print}{summary.crr}(x, digits = max(options()$digits - 4, 3), ...)
+}
+
+\arguments{
+ \item{object}{ An object of class crr (output from the crr function) }
+ \item{conf.int}{the level for a two-sided confidence interval on the
+ coeficients. Default is 0.95. }
+ \item{digits}{ In \code{summary.crr}, \code{digits} determines the number of
+ significant digits retained in the p-values. In
+ \code{print.summary.crr},
+ \code{digits} sets the values of the \code{digits} option for printing
+ the output.}
+ \item{\dots}{ Included for compatibility with the generic functions. Not
+ currently used. }
+ \item{x}{An object of class summary.crr (output from the summary
+ method for crr)}
+}
+\details{
+ The summary method calculates the standard errors, subdistribution
+ hazard ratios z-scores, p-values, and confidence intervals on the
+ hazard ratios. The print method prints a fairly standard format
+ tabular summary of the results.
+
+ The pseudo likelihood ratio test in the printed output is based on the
+ difference in the objective function at the global null and at the
+ final estimates. Since this objective function is not a true
+ likelihood, this test statistic is not asymptotically chi-square.
+}
+\value{\code{summary.crr} returns a list of class summary.crr, which
+ contains components
+ \item{call }{The call to crr}
+ \item{converged }{TRUE if the iterative algorithm converged}
+ \item{n }{The number of observations used in fitting the model}
+ \item{n.missing }{The number of observations removed by \code{crr}
+ from the input data due to missing values}
+ \item{loglik }{The value of the negative of the objective function
+ (the pseudo log likelihood at convergence}
+ \item{coef }{A matrix giving the estimated coefficients, hazard
+ ratios, standard errors, z-scores, and p-values }
+ \item{conf.int }{A matrix giving the estimated hazard ratios, inverse
+ hazard ratios and lower and upper confidence limits on the hazard ratios}
+ \item{logtest }{Twice the difference in log pseudo likelihood values}
+}
+\author{The summary and print.summary methods were provided by Luca Scrucca}
+
+\seealso{ \code{\link{crr}}}
+\examples{
+## see examples in the crr help file
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{survival}
diff --git a/man/timepoints.Rd b/man/timepoints.Rd
new file mode 100644
index 0000000..5546679
--- /dev/null
+++ b/man/timepoints.Rd
@@ -0,0 +1,34 @@
+\name{timepoints}
+\alias{timepoints}
+\title{
+Calculate Estimates at Specific Timepoints
+}
+\description{
+Find values at specified timepoints from curves specified as all
+corners of step functions.
+}
+\usage{
+timepoints(w, times)
+}
+\arguments{
+\item{w}{
+a list containing the estimates, with points for all corners of the
+step function. (Usually created by cuminc.)
+Each component in the list contains the estimate for a different group.
+Each component has components giving times, function estimates, and
+variances (see cuminc)
+}
+\item{times}{
+vector of times where estimates are needed
+}}
+\value{
+ A list with components
+ \item{$est}{a matrix of estimates of the subdistributions with a row for
+each component in \code{w} and a column for each time}
+ \item{$var}{a matrix giving the corresponding variances.}
+}
+\seealso{
+\code{\link{cuminc}}
+}
+\keyword{survival}
+% Converted by Sd2Rd version 1.10.
diff --git a/src/cincsub.f b/src/cincsub.f
new file mode 100644
index 0000000..91a2d5f
--- /dev/null
+++ b/src/cincsub.f
@@ -0,0 +1,104 @@
+c Copyright (C) 2000 Robert Gray
+c distributed under the terms of the GNU public license
+ subroutine cinc(y,ic,icc,n,x,f,v)
+c calculates estimates of cumulative incidence functions and their
+c variances
+c
+c i-n are integer, all others double precision
+c
+C DATA ARE ASSUMED TO BE SORTED (increasing order of y).
+c On input:
+c y is the failure times (assumed >= 0)
+c ic is 1 if the case has failed (from any cause), 0 otherwise
+c icc is 1 if the case has failed from the specific cause for which
+c the estimate is being calculated, 0 otherwise
+c n is the length of y (and ic and icc)
+c
+c x,f,v need to have length at least 2*nfc+2
+c where nfc is the # of unique failure times from the cause of
+c interest
+c On output:
+c x gives the times, f the estimates, and v the variances
+c This program was designed to produce estimates suitable for
+c drawing as a step function. x(0) is set to 0, and f(0) is also.
+c X(1) and x(2) are the value of the smallest time where icc=1
+c f(1) is 0, and f(2) is the value f jumps to at x(1)=x(2).
+c in this way all the corners of the step function are given in the
+c vector, and all the plotting routine has to do to make a plot is
+c to connect up the points.
+C THE LAST ELEMENT OF X AND F CONTNUE OUT TO THE largest follow-up time
+ implicit double precision (a-h,o-z)
+ dimension x(*),f(*),ic(n),icc(n),v(*),y(n)
+ fk=1
+ nf=0
+ v1=0
+ v2=0
+ v3=0
+ x(1)=0
+ f(1)=0
+ v(1)=0
+ lcnt=1
+ l=1
+ ll=1
+ rs=n
+ ty=y(1)
+ 10 l=l+1
+ if (l.gt.n) go to 60
+ if (y(l).eq.ty) go to 10
+ 60 l=l-1
+ nd1=0
+ nd2=0
+ do 15 i=ll,l
+ nd1=nd1+icc(i)
+ 15 nd2=nd2+ic(i)-icc(i)
+ nd=nd1+nd2
+ if (nd.eq.0) go to 40
+ fkn=fk*(rs-nd)/rs
+ if (nd1.gt.0) then
+ lcnt=lcnt+2
+ f(lcnt-1)=f(lcnt-2)
+ f(lcnt)=f(lcnt-1)+fk*nd1/rs
+ end if
+ if (nd2.le.0.or.fkn.le.0) go to 30
+ t5=1
+ if (nd2.gt.1) t5=1-(nd2-1.)/(rs-1.)
+ t6=fk*fk*t5*nd2/(rs*rs)
+ t3=1./fkn
+ t4=f(lcnt)/fkn
+ v1=v1+t4*t4*t6
+ v2=v2+t3*t4*t6
+ v3=v3+t3*t3*t6
+ 30 if (nd1.le.0) go to 35
+ t5=1
+ if (nd1.gt.1) t5=1-(nd1-1.)/(rs-1.)
+ t6=fk*fk*t5*nd1/(rs*rs)
+ t3=0
+ if (fkn.gt.0) t3=1/fkn
+ t4=1+t3*f(lcnt)
+ v1=v1+t4*t4*t6
+ v2=v2+t3*t4*t6
+ v3=v3+t3*t3*t6
+c changed 3-18-91: If largest follow-up time is a failure of type 1,
+c only v1 should be incremented (as above), but to get the right
+c variance still need to incorporate v2 and v3 in v(lcnt)
+c t2=0
+c if (fkn.gt.0) t2=f(lcnt)
+ t2=f(lcnt)
+ x(lcnt-1)=y(l)
+ x(lcnt)=y(l)
+ v(lcnt-1)=v(lcnt-2)
+ v(lcnt)=v1+t2*t2*v3-2*t2*v2
+ 35 fk=fkn
+ nf=nf+nd1
+ 40 rs=n-l
+ l=l+1
+ if (l.gt.n) go to 50
+ ll=l
+ ty=y(l)
+ go to 10
+ 50 lcnt=lcnt+1
+ x(lcnt)=y(n)
+ f(lcnt)=f(lcnt-1)
+ v(lcnt)=v(lcnt-1)
+ return
+ end
diff --git a/src/crr.f b/src/crr.f
new file mode 100644
index 0000000..f2f473e
--- /dev/null
+++ b/src/crr.f
@@ -0,0 +1,577 @@
+c Copyright (C) 2000 Robert Gray
+c distributed under the terms of the GNU public license
+ subroutine covt(j,k,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+c j = row index of case in x and x2
+c k = row index of cft in tf
+ double precision x(n,ncov),x2(n,ncov2),tf(ndf,ncov2),wk
+ double precision b(ncov+ncov2),xbt(ncov+ncov2)
+ integer ncov,ncov2,ndf,i,j,k,n
+ wk=0
+ if (ncov.gt.0) then
+ do 10 i=1,ncov
+ xbt(i)=x(j,i)
+ wk=wk+xbt(i)*b(i)
+ 10 continue
+ endif
+ if (ncov2.gt.0) then
+ do 11 i=1,ncov2
+ xbt(i+ncov)=x2(j,i)*tf(k,i)
+ wk=wk+xbt(i+ncov)*b(ncov+i)
+ 11 continue
+ endif
+ return
+ end
+
+ subroutine crrfsv(t2,ici,n,x,ncov,np,x2,ncov2,tf,ndf,wt,ncg,icg,b,
+ $ lik,s,v,xb,xbt,vt)
+c all data sorted with t2 in ascending order
+c t2, event/censoring time for each subject
+c ici, ici(i)=1 if t2(i) is a type 1 failure time (& =2 for other failures)
+c & =0 for censored
+c x(n,ncov) ph covariates, x2(n,ncov2) covs multiplied by functions of time
+c tf(i,j) is the value of the time function which multiplies the jth col of
+c x2, at the ith distinct type 1 failure time (ascending order);
+c ndf is the number of distinct type 1 failures=# rows in tf
+c wt are km est of censoring dists at times in t2- for each group formed by
+c distinct values of icg.
+c censoring groups in icg must be coded 1,2,...,ncg
+c output in lik,s,v
+c xb,xbt,vt are temporary storage
+ double precision t2(n),b(np),s(np),v(np,np),lik,tf(ndf,ncov2)
+ double precision xb(np),xb1,vt(np,np),wt(ncg,n),cft,twf
+ double precision x(n,ncov),x2(n,ncov2),wk,xb1o,twt,xbt(np)
+ integer ici(n),icg(n),n,np,i,j,iuc,ncov,ncov2,ncg,ndf,ldf,k,itmp
+ lik=0.d0
+ do 1 i=1,np
+ s(i)=0
+ do 2 j=i,np
+ v(i,j)=0
+ 2 continue
+ 1 continue
+c
+c to see how this works, use the simple nested loop approach
+c
+ iuc=n
+ ldf=ndf+1
+c find next failure time
+ 98 itmp=iuc
+ do 10 i=iuc,1,-1
+ itmp=i
+ if (ici(i).eq.1) then
+ cft=t2(i)
+ go to 11
+ endif
+ 10 continue
+c no more failures
+ return
+ 11 iuc=itmp
+ ldf=ldf-1
+ twf=0
+ do 13 i=iuc,1,-1
+ if (t2(i).lt.cft) go to 14
+ itmp=i
+ if (ici(i).eq.1) then
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twf=twf+1
+c using a minimization routine, so neg of objective fcn, scores, etc
+ lik=lik-wk
+ do 19 j=1,np
+ s(j)=s(j)-xbt(j)
+ 19 continue
+ endif
+ 13 continue
+c calculate sums over risk set
+ 14 iuc=itmp
+ xb1=0
+ xb1o=xb1
+ do 3 i=1,np
+ xb(i)=0
+ do 4 j=i,np
+ vt(i,j)=0
+ 4 continue
+ 3 continue
+ do 15 i=1,n
+ if (t2(i).lt.cft) then
+ if (ici(i).le.1) go to 15
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(i),iuc)/wt(icg(i),i)
+ else
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ endif
+ xb1=xb1+twt
+ do 17 j=1,np
+ xb(j)=xb(j)+twt*xbt(j)
+ xbt(j)=xbt(j)-xb(j)/xb1
+ 17 continue
+ if (xb1o.gt.0) then
+ twt=xb1*twt/xb1o
+ do 51 k=1,np
+ do 52 j=k,np
+ vt(k,j)=vt(k,j)+twt*xbt(k)*xbt(j)
+ 52 continue
+ 51 continue
+ endif
+c all dsyr('U',np,xb1*twt/xb1o,xbt,1,vt,np)
+ xb1o=xb1
+ 15 continue
+ lik=lik+twf*log(xb1)
+ twt=twf/xb1
+ do 61 i=1,np
+ s(i)=s(i)+twt*xb(i)
+ do 62 j=i,np
+ v(i,j)=v(i,j)+twt*vt(i,j)
+ v(j,i)=v(i,j)
+ 62 continue
+ 61 continue
+c call daxpy(np,-twf/xb1,xb,1,s,1)
+c call daxpy(np*np,twf/xb1,vt,1,v,1)
+ iuc=iuc-1
+ if (iuc.gt.0) go to 98
+ return
+ end
+
+ subroutine crrf(t2,ici,n,x,ncov,np,x2,ncov2,tf,ndf,wt,ncg,icg,b,
+ $ lik,xbt)
+c all data sorted with t2 in ascending order
+c t2, stop time for each subject
+c ici, ici(i)=1 if t2(i) is a type 1 failure time (& =2 for other failures)
+c x(nrx,np) covariates
+c wt are km est of censoring dists at times in t2-
+c output in lik,s,v
+c icrs,wk,xb,xbt,vt are temporary storage
+ double precision t2(n),b(np),lik,tf(ndf,ncov2)
+ double precision xb1,wt(ncg,n),cft,twf,xbt(np)
+ double precision x(n,ncov),x2(n,ncov2),wk,twt
+ integer ici(n),icg(n),n,np,i,iuc,ncov,ncov2,ncg,ndf,ldf,itmp
+ lik=0.d0
+ iuc=n
+ ldf=ndf+1
+c find next failure time
+ 98 itmp=iuc
+ do 10 i=iuc,1,-1
+ itmp=i
+ if (ici(i).eq.1) then
+ cft=t2(i)
+ go to 11
+ endif
+ 10 continue
+c no more failures
+ return
+ 11 iuc=itmp
+ ldf=ldf-1
+ twf=0
+ do 13 i=iuc,1,-1
+ if (t2(i).lt.cft) go to 14
+ itmp=i
+ if (ici(i).eq.1) then
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twf=twf+1
+ lik=lik-wk
+ endif
+ 13 continue
+c calculate sums over risk set
+ 14 iuc=itmp
+ xb1=0
+ do 15 i=1,n
+ if (t2(i).lt.cft) then
+ if (ici(i).le.1) go to 15
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(i),iuc)/wt(icg(i),i)
+ else
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ endif
+ xb1=xb1+twt
+ 15 continue
+ lik=lik+twf*log(xb1)
+ iuc=iuc-1
+ if (iuc.gt.0) go to 98
+ return
+ end
+
+ subroutine crrvv(t2,ici,n,x,ncov,np,x2,ncov2,tf,ndf,wt,ncg,icg,b,
+ $ v,v2,vt,xb,xbt,qu,st,ss2,icrsk,ss3,ss4)
+c all data sorted with t2 in ascending order
+c t2, event/censoring time for each subject
+c ici, ici(i)=1 if t2(i) is a type 1 failure time (& =2 for other failures)
+c & =0 for censored
+c x(n,ncov) ph covariates, x2(n,ncov2) covs multiplied by functions of time
+c tf(i,j) is the value of the time function which multiplies the jth col of
+c x2, at the ith distinct type 1 failure time (ascending order);
+c ndf is the number of distinct type 1 failures=# rows in tf
+c wt are km est of censoring dists at times in t2- for each group formed by
+c distinct values of icg.
+c censoring groups in icg must be coded 1,2,...,ncg
+c output in v,v2
+c xb,xbt,vt,ss,st,qu,icrsk are temporary storage
+ double precision t2(n),b(np),v(np,np),tf(ndf,ncov2)
+ double precision xb(n,0:np),vt(np,np),wt(ncg,n),cft
+ double precision x(n,ncov),x2(n,ncov2),wk,twt,xbt(np),ss2(np,ncg)
+ double precision st(np,2),v2(np,np),cft2,qu(np,ncg)
+ double precision ss3(np,ncg),ss4(ncg)
+ integer ici(n),icg(n),n,np,i,j,ncov,ncov2,ncg,ndf,ldf,k
+ integer lc,icrsk(ncg),j1,j2,ldf2,iflg
+ do 3 i=1,ncg
+ icrsk(i)=0
+ 3 continue
+ do 1 i=1,np
+ do 123 j=1,ncg
+ ss2(i,j)=0
+ 123 continue
+ do 2 j=i,np
+ v(i,j)=0
+ v2(i,j)=0
+ 2 continue
+ 1 continue
+ do 5 i=1,n
+ icrsk(icg(i))=icrsk(icg(i))+1
+ do 4 j=0,np
+ xb(i,j)=0
+ 4 continue
+ 5 continue
+ ldf=0
+ cft=min(-1.d0,t2(1)*(1-1.d-5))
+ do 6 i=1,n
+ if (ici(i).ne.1) go to 6
+ if (t2(i).gt.cft) then
+ cft=t2(i)
+ ldf=ldf+1
+ endif
+ do 7 j=1,n
+ if (t2(j).lt.t2(i)) then
+ if (ici(j).le.1) go to 7
+ call covt(j,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(j),i)/wt(icg(j),j)
+ else
+ call covt(j,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ endif
+ xb(i,0)=xb(i,0)+twt
+ do 8 k=1,np
+ xb(i,k)=xb(i,k)+twt*xbt(k)
+ 8 continue
+ 7 continue
+ 6 continue
+c
+ lc=1
+ ldf2=0
+ cft2=min(-1.d0,t2(1)*(1-1.d-5))
+ do 10 i=1,n
+ do 11 k=1,np
+ st(k,1)=0
+ 11 continue
+ ldf=0
+ cft=min(-1.d0,t2(1)*(1-1.d-5))
+ do 15 j=1,n
+ if (ici(j).ne.1) go to 15
+ if (t2(j).gt.cft) then
+ cft=t2(j)
+ ldf=ldf+1
+ endif
+ if (t2(j).le.t2(i)) then
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ else if (t2(i).lt.t2(j).and.ici(i).gt.1) then
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(i),j)/wt(icg(i),i)
+ else
+ go to 15
+ endif
+c d lambda hat portion of eta_i
+ do 16 k=1,np
+ st(k,1)=st(k,1)-(xbt(k)-xb(j,k)/xb(j,0))*twt/xb(j,0)
+ 16 continue
+ 15 continue
+ if (ici(i).eq.1) then
+ if (t2(i).gt.cft2) then
+ cft2=t2(i)
+ ldf2=ldf2+1
+ endif
+ call covt(i,ldf2,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+c d N_i portion of eta_i
+ do 17 k=1,np
+ st(k,1)=st(k,1)+xbt(k)-xb(i,k)/xb(i,0)
+ 17 continue
+c second derivatives
+ do 201 j1=1,np
+ do 202 j2=j1,np
+ vt(j1,j2)=0
+ 202 continue
+ 201 continue
+ do 19 j=1,n
+ if (t2(j).lt.t2(i)) then
+ if (ici(j).le.1) go to 19
+ call covt(j,ldf2,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(j),i)/wt(icg(j),j)
+ else
+ call covt(j,ldf2,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ endif
+ do 18 k=1,np
+ xbt(k)=xbt(k)-xb(i,k)/xb(i,0)
+ 18 continue
+ do 51 k=1,np
+ do 52 j2=k,np
+ vt(k,j2)=vt(k,j2)+twt*xbt(k)*xbt(j2)
+ 52 continue
+ 51 continue
+ 19 continue
+ do 251 j1=1,np
+ do 252 j2=j1,np
+ v(j1,j2)=v(j1,j2)+vt(j1,j2)/xb(i,0)
+ 252 continue
+ 251 continue
+ endif
+c call dblepr('st1',3,st(1,1),np)
+ if (i.eq.1) then
+ iflg=1
+ else
+ if (t2(i).gt.t2(i-1)) then
+ iflg=1
+ else
+ iflg=0
+ endif
+ endif
+C if (i.eq.1.or.t2(i).gt.t2(i-1)) then
+ if (iflg.eq.1) then
+ do 40 j=i,n
+ if (t2(j).gt.t2(i)) go to 39
+ if (ici(j).eq.0) go to 38
+ 40 continue
+c calculate qu (q(u))
+ 38 ldf=ldf2
+ cft=cft2
+ do 280 k2=1,ncg
+ do 281 k=1,np
+ qu(k,k2)=0
+ 281 continue
+ 280 continue
+c j1 indexes inner integral in q (t2(i) is lower limit of integration)
+c dN_j2 part is 0, because integrand includes I(s>t2(j2)), where s is var
+c of integration, so the sum over j1 is calculating the d lambda_1 hat part
+ do 41 j1=lc,n
+ if (ici(j1).ne.1) go to 41
+ if (t2(j1).gt.cft) then
+ cft=t2(j1)
+ ldf=ldf+1
+ endif
+ do 555 k3=1,ncg
+ ss4(k3)=0
+ do 255 k=1,np
+ ss3(k,k3)=0
+ 255 continue
+ 555 continue
+c j2 indexes outer sum in q. because of I(s>t2(j2)) and
+c I(t2(j2)>=s)+I(t2(j2)<s,eps_j2=2), only type 2 failures contribute
+ do 541 j2=1,n
+ if (t2(j2).ge.t2(i)) go to 542
+ if (ici(j2).le.1) go to 541
+ call covt(j2,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(j2),j1)/wt(icg(j2),j2)
+ ss4(icg(j2))=ss4(icg(j2))+twt
+ do 256 k=1,np
+ ss3(k,icg(j2))=ss3(k,icg(j2))+xbt(k)*twt
+ 256 continue
+ 541 continue
+c qu is q(t2(i))
+ 542 do 42 k=1,np
+ qu(k,icg(j1))=qu(k,icg(j1))+(ss3(k,icg(j1))-xb(j1,k)*
+ $ ss4(icg(j1))/xb(j1,0))/xb(j1,0)
+ 42 continue
+ 41 continue
+c call dblepr('qu',2,qu(1),1)
+c update ss2 (note the sign, this is already negative)--ss2 is the
+c d lambda c hat portion of psi_i. This is the addition
+c to the d lambda c hat portion at the ith censoring time.
+ do 43 j=i,n
+ if (t2(j).gt.t2(i)) go to 39
+ if (ici(j).eq.0) then
+ do 282 k=1,np
+ ss2(k,icg(j))=ss2(k,icg(j))-qu(k,icg(j))/
+ $ icrsk(icg(j))**2
+ 282 continue
+ endif
+ 43 continue
+ endif
+c st(,2) is psi
+ 39 do 243 k=1,np
+ st(k,2)=ss2(k,icg(i))
+ 243 continue
+c dNc portion of psi
+ if (ici(i).eq.0) then
+ do 45 k=1,np
+ st(k,2)=st(k,2)+qu(k,icg(i))/icrsk(icg(i))
+ 45 continue
+ endif
+ do 271 j1=np,1,-1
+ st(j1,1)=st(j1,1)+st(j1,2)
+ do 272 j2=j1,np
+ v2(j1,j2)=v2(j1,j2)+st(j1,1)*st(j2,1)
+ 272 continue
+ 271 continue
+c call dblepr('st1',3,st(1,1),np)
+ if (i.lt.n) then
+ if (t2(i+1).gt.t2(i)) then
+ do 249 j=lc,i
+ icrsk(icg(j))=icrsk(icg(j))-1
+ 249 continue
+ lc=i+1
+ endif
+ endif
+ 10 continue
+c call intpr('icrsk',5,icrsk,2)
+ do 333 j1=1,np-1
+ do 334 j2=j1+1,np
+ v(j2,j1)=v(j1,j2)
+ v2(j2,j1)=v2(j1,j2)
+ 334 continue
+ 333 continue
+ return
+ end
+
+ subroutine crrsr(t2,ici,n,x,ncov,np,x2,ncov2,tf,ndf,wt,ncg,icg,b,
+ $ res,xb,xbt)
+c compute contribution to scores at each type 1 failure time
+c all data sorted with t2 in ascending order
+c t2, event/censoring time for each subject
+c ici, ici(i)=1 if t2(i) is a type 1 failure time (& =2 for other failures)
+c & =0 for censored
+c x(n,ncov) ph covariates, x2(n,ncov2) covs multiplied by functions of time
+c tf(i,j) is the value of the time function which multiplies the jth col of
+c x2, at the ith distinct type 1 failure time (ascending order);
+c ndf is the number of distinct type 1 failures=# rows in tf
+c wt are km est of censoring dists at times in t2- for each group formed by
+c distinct values of icg.
+c score residuals in res
+c xb,xbt are temporary storage
+ double precision t2(n),b(np),tf(ndf,ncov2)
+ double precision xb(np),wt(ncg,n),cft,res(np,ndf)
+ double precision x(n,ncov),x2(n,ncov2),wk,twt,twf,xbt(np),xb1
+ integer ici(n),icg(n),n,np,i,j,ncov,ncov2,ncg,ndf,ldf,iuc,itmp
+ do 1 i=1,ndf
+ do 2 j=1,np
+ res(j,i)=0
+ 2 continue
+ 1 continue
+ iuc=n
+ ldf=ndf+1
+c find next failure time
+ 98 itmp=iuc
+ do 10 i=iuc,1,-1
+ itmp=i
+ if (ici(i).eq.1) then
+ cft=t2(i)
+ go to 11
+ endif
+ 10 continue
+c no more failures
+ return
+ 11 iuc=itmp
+ ldf=ldf-1
+ twf=0
+ do 13 i=iuc,1,-1
+ if (t2(i).lt.cft) go to 14
+ itmp=i
+ if (ici(i).eq.1) then
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twf=twf+1
+ do 241 j=1,np
+ res(j,ldf)=res(j,ldf)+xbt(j)
+ 241 continue
+ endif
+ 13 continue
+c calculate sums over risk set
+ 14 iuc=itmp
+ xb1=0
+ do 3 i=1,np
+ xb(i)=0
+ 3 continue
+ do 15 i=1,n
+ if (t2(i).lt.cft) then
+ if (ici(i).le.1) go to 15
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(i),iuc)/wt(icg(i),i)
+ else
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ endif
+ xb1=xb1+twt
+ do 17 j=1,np
+ xb(j)=xb(j)+twt*xbt(j)
+ 17 continue
+ 15 continue
+ twt=-twf/xb1
+ do 61 i=1,np
+ res(i,ldf)=res(i,ldf)+twt*xb(i)
+ 61 continue
+ iuc=iuc-1
+c call intpr('iuc',3,iuc,1)
+ if (iuc.gt.0) go to 98
+ return
+ end
+
+ subroutine crrfit(t2,ici,n,x,ncov,np,x2,ncov2,tf,ndf,wt,ncg,icg,b,
+ $ res,xbt)
+c computes jumps in estimate of cumulative underlying hazard
+c at each type 1 failure time
+c all data sorted with t2 in ascending order
+c t2, event/censoring time for each subject
+c ici, ici(i)=1 if t2(i) is a type 1 failure time (& =2 for other failures)
+c & =0 for censored
+c x(n,ncov) ph covariates, x2(n,ncov2) covs multiplied by functions of time
+c tf(i,j) is the value of the time function which multiplies the jth col of
+c x2, at the ith distinct type 1 failure time (ascending order);
+c ndf is the number of distinct type 1 failures=# rows in tf
+c wt are km est of censoring dists at times in t2- for each group formed by
+c distinct values of icg.
+c res(i) is the jump in the estimated underlying cum pseudo haz for type 1
+c failures at the ith distinct (ascending) type 1 failure time
+c xbt is temporary storage
+ double precision t2(n),b(np),tf(ndf,ncov2)
+ double precision wt(ncg,n),cft,res(ndf)
+ double precision x(n,ncov),x2(n,ncov2),wk,twt,twf,xbt(np),xb1
+ integer ici(n),icg(n),n,np,i,ncov,ncov2,ncg,ndf,ldf,iuc,itmp
+ do 1 i=1,ndf
+ res(i)=0
+ 1 continue
+ iuc=1
+ ldf=0
+c find next failure time
+ 98 itmp=iuc
+ do 10 i=iuc,n
+ itmp=i
+ if (ici(i).eq.1) then
+ cft=t2(i)
+ go to 11
+ endif
+ 10 continue
+c no more failures
+ return
+ 11 iuc=itmp
+ ldf=ldf+1
+ twf=0
+ do 13 i=iuc,n
+ if (t2(i).gt.cft) go to 14
+ itmp=i
+ if (ici(i).eq.1) twf=twf+1
+ 13 continue
+c calculate sums over risk set
+ 14 iuc=itmp
+ xb1=0
+ do 15 i=1,n
+ if (t2(i).lt.cft) then
+ if (ici(i).le.1) go to 15
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)*wt(icg(i),iuc)/wt(icg(i),i)
+ else
+ call covt(i,ldf,ncov,x,n,ncov2,x2,tf,ndf,b,wk,xbt)
+ twt=exp(wk)
+ endif
+ xb1=xb1+twt
+ 15 continue
+ res(ldf)=res(ldf)+twf/xb1
+ iuc=iuc+1
+ if (iuc.le.n) go to 98
+ return
+ end
diff --git a/src/crstm.f b/src/crstm.f
new file mode 100644
index 0000000..72a45f9
--- /dev/null
+++ b/src/crstm.f
@@ -0,0 +1,244 @@
+c Copyright (C) 2000 Robert Gray
+c distributed under the terms of the GNU public license
+ subroutine crstm(y,m,ig,ist,no,rho,nst,ng,s,vs,ys,
+ &ms,igs,v,st,vt,wk,iwk)
+c
+c subroutine to calculate score and variance matrix for comparing
+c cumulative incidence curves for a specific cause among groups.
+c test statistic given by s' inv(vs) s, dist approx chi-square(ng-1)
+c
+c everything starting with i-n is integer, all others double precision
+c
+c On input:
+c y is the failure times (sorted in increasing order)
+c m is coded 0 if censored, 1 if failed from the cause of interest
+c 2 if failed from some other cause.
+c ig denotes group membership, must be coded 1,2,...,ng (ie
+c consecutive integers from 1 to ng), where ng is the number of groups
+c ist denotes strata membership, must be coded 1,2,...,nst, where nst
+c is the number of strata (code all 1's if only 1 strata)
+c no is the # of observations (length of y,m,ig,ist)
+c rho is the power used in the weight function in the test statistic
+c nst and ng are the # strata and # groups
+c Length of ys, ms, and igs must be at least as long as the size of the
+c largest strata
+c Length of s and st must be at least ng-1
+c Length of v and vt must be at least ng*(ng-1)/2
+c Length of vs must be at least (ng-1)^2
+c Length of wk must be at least ng*(4+3*ng)
+c Length of iwk must be at least 4*ng
+c
+c On output:
+c s gives the scores for the first ng-1 groups, and
+c vs the estimated variance covariance matrix of these scores
+c
+ implicit double precision (a-h,o-z)
+ dimension y(no),m(no),ig(no),ist(no),ys(no),ms(no)
+ dimension wk(ng*(4+3*ng)),iwk(4*ng)
+ dimension igs(no),s(ng-1),v(ng*(ng-1)/2),st(ng-1),vt(ng*(ng-1)/2)
+ dimension vs(ng-1,ng-1)
+ ng1=ng-1
+ ng2=ng*ng1/2
+ l=0
+ do 12 i=1,ng1
+ s(i)=0
+ do 14 j=1,i
+ l=l+1
+ v(l)=0
+ 14 continue
+ 12 continue
+ do 20 ks=1,nst
+ n=0
+ do 21 i=1,no
+ if (ist(i).ne.ks) go to 21
+ n=n+1
+ ys(n)=y(i)
+ ms(n)=m(i)
+ igs(n)=ig(i)
+ 21 continue
+ ng3=4*ng+1
+ ng4=ng*ng
+ call crst(ys(1),ms(1),igs(1),n,ng,rho,st,vt,ng1,ng2,wk(1),
+ &wk(ng+1),wk(2*ng+1),wk(3*ng+1),wk(ng3),wk(ng3+ng4),wk(ng3+2*ng4),
+ &wk(ng3+2*ng4+ng),iwk(1),iwk(ng+1))
+ l=0
+ do 23 i=1,ng1
+ s(i)=s(i)+st(i)
+ do 24 j=1,i
+ l=l+1
+ v(l)=v(l)+vt(l)
+ 24 continue
+ 23 continue
+ 20 continue
+ l=0
+ do 31 i=1,ng1
+ do 332 j=1,i
+ l=l+1
+ vs(i,j)=v(l)
+ vs(j,i)=vs(i,j)
+ 332 continue
+ 31 continue
+ return
+ end
+
+ subroutine crst(y,m,ig,n,ng,rho,s,v,ng1,nv,f1m,f1,skmm,
+ &skm,c,a,v3,v2,rs,d)
+ implicit double precision (a-h,o-z)
+ dimension y(n),s(ng1),f1m(ng),f1(ng),skmm(ng),skm(ng)
+ dimension c(ng,ng),a(ng,ng),v(nv),v3(ng)
+ dimension v2(ng1,ng)
+ integer m(n),ig(n),rs(ng),d(0:2,ng)
+c rs(j) will be the risk set size in group j at the current failure time
+c (initially the sample size in each group)
+ do 15 i=1,ng
+ 15 rs(i)=0
+ do 16 i=1,n
+ j=ig(i)
+ 16 rs(j)=rs(j)+1
+ l=0
+ do 11 i=1,ng1
+ s(i)=0
+ do 12 j=1,i
+ l=l+1
+ 12 v(l)=0
+ 11 continue
+ do 14 i=1,ng
+ f1m(i)=0
+ f1(i)=0
+ skmm(i)=1
+ skm(i)=1
+ v3(i)=0
+ do 9 j=1,ng1
+ 9 v2(j,i)=0
+ do 8 j=1,ng
+ 8 c(i,j)=0
+ 14 continue
+ fm=0
+ f=0
+c begin looping over unique times:
+ ll=1
+ lu=ll
+ 50 lu=lu+1
+ if (lu.gt.n) go to 55
+ if (y(lu).gt.y(ll)) go to 55
+ go to 50
+ 55 lu=lu-1
+ nd1=0
+ nd2=0
+c d will contain the # in each group censored, failed from
+c cause 1, and failing from cause 2, at this time
+ do 56 i=1,ng
+ d(0,i)=0
+ d(1,i)=0
+ 56 d(2,i)=0
+ do 57 i=ll,lu
+ j=ig(i)
+ k=m(i)
+ 57 d(k,j)=d(k,j)+1
+ do 58 i=1,ng
+ nd1=nd1+d(1,i)
+ nd2=nd2+d(2,i)
+ 58 continue
+ if (nd1.eq.0.and.nd2.eq.0) go to 90
+ tr=0
+ tq=0
+ do 60 i=1,ng
+ if (rs(i).le.0) go to 60
+ td=d(1,i)+d(2,i)
+c skmm is left continuous, and skm right continuous, km est.
+ skm(i)=skmm(i)*(rs(i)-td)/rs(i)
+c f1m is left continuous, and f1 right continuous, cuminc est.
+ f1(i)=f1m(i)+(skmm(i)*d(1,i))/rs(i)
+c in notation of the paper, tr is \sum_r\hat{h}_r, and tq is \sum_r R_r
+ tr=tr+rs(i)/skmm(i)
+ tq=tq+rs(i)*(1-f1m(i))/skmm(i)
+ 60 continue
+ f=fm+nd1/tr
+ fb=(1-fm)**rho
+ do 66 i=1,ng
+ do 166 j=i,ng
+ 166 a(i,j)=0
+ if (rs(i).le.0) go to 66
+ t1=rs(i)/skmm(i)
+ a(i,i)=fb*t1*(1-t1/tr)
+ c(i,i)=c(i,i)+a(i,i)*nd1/(tr*(1-fm))
+ k=i+1
+ if (k.gt.ng) go to 66
+ do 67 j=k,ng
+ if (rs(j).le.0) go to 67
+ a(i,j)=-fb*t1*rs(j)/(skmm(j)*tr)
+ c(i,j)=c(i,j)+a(i,j)*nd1/(tr*(1-fm))
+ 67 continue
+ 66 continue
+ do 68 i=2,ng
+ k=i-1
+ do 69 j=1,k
+ a(i,j)=a(j,i)
+ 69 c(i,j)=c(j,i)
+ 68 continue
+ do 74 i=1,ng1
+ if (rs(i).le.0) go to 74
+ s(i)=s(i)+fb*(d(1,i)-nd1*rs(i)*(1-f1m(i))/(skmm(i)*tq))
+ 74 continue
+ if (nd1.le.0) go to 77
+ do 72 k=1,ng
+ if (rs(k).le.0) go to 72
+ t4=1
+ if (skm(k).gt.0) t4=1-(1-f)/skm(k)
+ t5=1
+ if (nd1.gt.1) t5=1-(nd1-1)/(tr*skmm(k)-1)
+ t3=t5*skmm(k)*nd1/(tr*rs(k))
+ v3(k)=v3(k)+t4*t4*t3
+ do 70 i=1,ng1
+ t1=a(i,k)-t4*c(i,k)
+ v2(i,k)=v2(i,k)+t1*t4*t3
+ do 71 j=1,i
+ l=i*(i-1)/2+j
+ t2=a(j,k)-t4*c(j,k)
+ v(l)=v(l)+t1*t2*t3
+ 71 continue
+ 70 continue
+ 72 continue
+ 77 if (nd2.eq.0) go to 90
+ do 82 k=1,ng
+ if (skm(k).le.0.or.d(2,k).le.0) go to 82
+ t4=(1-f)/skm(k)
+ t5=1
+c following line changed 3-24-04 - had been performed as integer
+ if (d(2,k).gt.1) t5=1-(d(2,k)-1.d0)/(rs(k)-1.d0)
+ t3=t5*((skmm(k)**2)*d(2,k))/(rs(k)**2)
+ v3(k)=v3(k)+t4*t4*t3
+ do 80 i=1,ng1
+ t1=t4*c(i,k)
+ v2(i,k)=v2(i,k)-t1*t4*t3
+ do 81 j=1,i
+ l=i*(i-1)/2+j
+ t2=t4*c(j,k)
+ v(l)=v(l)+t1*t2*t3
+ 81 continue
+ 80 continue
+ 82 continue
+ 90 if (lu.ge.n) go to 30
+ do 91 i=ll,lu
+ j=ig(i)
+ 91 rs(j)=rs(j)-1
+ fm=f
+ do 92 i=1,ng
+ f1m(i)=f1(i)
+ 92 skmm(i)=skm(i)
+ ll=lu+1
+ lu=ll
+ go to 50
+ 30 l=0
+ do 36 i=1,ng1
+ do 37 j=1,i
+ l=l+1
+ do 38 k=1,ng
+ v(l)=v(l)+c(i,k)*c(j,k)*v3(k)
+ v(l)=v(l)+c(i,k)*v2(j,k)
+ v(l)=v(l)+c(j,k)*v2(i,k)
+ 38 continue
+ 37 continue
+ 36 continue
+ return
+ end
diff --git a/src/tpoi.f b/src/tpoi.f
new file mode 100644
index 0000000..47243b8
--- /dev/null
+++ b/src/tpoi.f
@@ -0,0 +1,34 @@
+c Copyright (C) 2000 Robert Gray
+c Distributed under the terms of the GNU public license
+ subroutine tpoi(x,n,ind,tp5,ntp)
+ double precision x(n),tp5(ntp)
+ integer n,ntp,ind(ntp),l,k,i
+ l=ntp
+ do 10 i=ntp,1,-1
+ if (x(n).ge.tp5(i)) go to 9
+ ind(l)=0
+ l=l-1
+ 10 continue
+ 9 if (l.le.0) return
+ if (x(n).eq.tp5(l)) then
+ ind(l)=n
+ l=l-1
+ endif
+c assuming unique values sorted in ascending order
+ k=n-1
+ 11 if (l.le.0) return
+ do 20 i=k,1,-1
+ if (x(k).le.tp5(l)) then
+ ind(l)=k+1
+ l=l-1
+ go to 11
+ else
+ k=k-1
+ endif
+ 20 continue
+c error in the following loop corrected 9-28-04
+ do 30 i=1,l
+ ind(i)=0
+ 30 continue
+ return
+ end
diff --git a/tests/Rplots.ps b/tests/Rplots.ps
new file mode 100644
index 0000000..0466870
--- /dev/null
+++ b/tests/Rplots.ps
@@ -0,0 +1,1340 @@
+%!PS-Adobe-3.0
+%%DocumentNeededResources: font Helvetica
+%%+ font Helvetica-Bold
+%%+ font Helvetica-Oblique
+%%+ font Helvetica-BoldOblique
+%%+ font Symbol
+%%DocumentMedia: a4 595 841 0 () ()
+%%Title: R Graphics Output
+%%Creator: R Software
+%%Pages: (atend)
+%%Orientation: Landscape
+%%BoundingBox: 18 18 577 824
+%%EndComments
+%%BeginProlog
+/bp { gs 595.00 0 translate 90 rotate gs } def
+% begin .ps.prolog
+/gs { gsave } def
+/gr { grestore } def
+/ep { showpage gr gr } def
+/m { moveto } def
+/l { rlineto } def
+/np { newpath } def
+/cp { closepath } def
+/f { fill } def
+/o { stroke } def
+/c { newpath 0 360 arc } def
+/r { 4 2 roll moveto 1 copy 3 -1 roll exch 0 exch rlineto 0 rlineto -1 mul 0 exch rlineto closepath } def
+/p1 { stroke } def
+/p2 { gsave bg fill grestore newpath } def
+/p3 { gsave bg fill grestore stroke } def
+/t { 6 -2 roll moveto gsave rotate
+ ps mul neg 0 2 1 roll rmoveto
+ 1 index stringwidth pop
+ mul neg 0 rmoveto show grestore } def
+/cl { grestore gsave newpath 3 index 3 index moveto 1 index
+ 4 -1 roll lineto exch 1 index lineto lineto
+ closepath clip newpath } def
+/rgb { setrgbcolor } def
+/s { scalefont setfont } def
+% end .ps.prolog
+%%IncludeResource: font Helvetica
+/Helvetica findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding ISOLatin1Encoding def
+ currentdict
+ end
+/Font1 exch definefont pop
+%%IncludeResource: font Helvetica-Bold
+/Helvetica-Bold findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding ISOLatin1Encoding def
+ currentdict
+ end
+/Font2 exch definefont pop
+%%IncludeResource: font Helvetica-Oblique
+/Helvetica-Oblique findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding ISOLatin1Encoding def
+ currentdict
+ end
+/Font3 exch definefont pop
+%%IncludeResource: font Helvetica-BoldOblique
+/Helvetica-BoldOblique findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ /Encoding ISOLatin1Encoding def
+ currentdict
+ end
+/Font4 exch definefont pop
+%%IncludeResource: font Symbol
+/Symbol findfont
+dup length dict begin
+ {1 index /FID ne {def} {pop pop} ifelse} forall
+ currentdict
+ end
+/Font5 exch definefont pop
+%%EndProlog
+%%Page: 1 1
+bp
+18.00 18.00 823.89 577.28 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 91.44 m
+632.09 0 l
+o
+np
+103.58 91.44 m
+0 -7.20 l
+o
+np
+230.00 91.44 m
+0 -7.20 l
+o
+np
+356.42 91.44 m
+0 -7.20 l
+o
+np
+482.84 91.44 m
+0 -7.20 l
+o
+np
+609.26 91.44 m
+0 -7.20 l
+o
+np
+735.67 91.44 m
+0 -7.20 l
+o
+/ps 12 def /Font1 findfont 12 s
+103.58 65.52 (0) .5 0 0 t
+230.00 65.52 (1) .5 0 0 t
+356.42 65.52 (2) .5 0 0 t
+482.84 65.52 (3) .5 0 0 t
+609.26 65.52 (4) .5 0 0 t
+735.67 65.52 (5) .5 0 0 t
+np
+77.04 107.25 m
+0 395.18 l
+o
+np
+77.04 107.25 m
+-7.20 0 l
+o
+np
+77.04 186.28 m
+-7.20 0 l
+o
+np
+77.04 265.32 m
+-7.20 0 l
+o
+np
+77.04 344.36 m
+-7.20 0 l
+o
+np
+77.04 423.39 m
+-7.20 0 l
+o
+np
+77.04 502.43 m
+-7.20 0 l
+o
+59.76 107.25 (0.0) .5 0 90 t
+59.76 186.28 (0.2) .5 0 90 t
+59.76 265.32 (0.4) .5 0 90 t
+59.76 344.36 (0.6) .5 0 90 t
+59.76 423.39 (0.8) .5 0 90 t
+59.76 502.43 (1.0) .5 0 90 t
+np
+77.04 518.24 m
+0 -426.80 l
+716.61 0 l
+o
+18.00 18.00 823.89 577.28 cl
+/ps 14 def /Font2 findfont 14 s
+0 setgray
+435.34 542.73 ( ) .5 0 0 t
+/ps 12 def /Font1 findfont 12 s
+435.34 36.72 (Years) .5 0 0 t
+30.96 304.84 (Probability) .5 0 90 t
+77.04 91.44 793.65 518.24 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+114.38 488.03 m
+21.60 0 l
+o
+0.75 setlinewidth
+[ 3.00 5.00] 0 setdash
+np
+114.38 473.63 m
+21.60 0 l
+o
+0.75 setlinewidth
+[ 0.00 4.00] 0 setdash
+np
+114.38 459.23 m
+21.60 0 l
+o
+0.75 setlinewidth
+[ 0.00 4.00 3.00 4.00] 0 setdash
+np
+114.38 444.83 m
+21.60 0 l
+o
+0.75 setlinewidth
+[ 6.00 4.00] 0 setdash
+np
+114.38 430.43 m
+21.60 0 l
+o
+0.75 setlinewidth
+[ 1.00 3.00 5.00 3.00] 0 setdash
+np
+114.38 416.03 m
+21.60 0 l
+o
+/ps 12 def /Font1 findfont 12 s
+146.78 483.72 (a 1) 0 0 0 t
+146.78 469.32 (b 1) 0 0 0 t
+146.78 454.92 (c 1) 0 0 0 t
+146.78 440.52 (a 2) 0 0 0 t
+146.78 426.12 (b 2) 0 0 0 t
+146.78 411.72 (c 2) 0 0 0 t
+0.75 setlinewidth
+[] 0 setdash
+np
+103.58 107.25 m
+17.20 0 l
+0 11.97 l
+22.03 0 l
+0 12.86 l
+26.56 0 l
+0 12.87 l
+30.66 0 l
+0 14.12 l
+30.46 0 l
+0 14.12 l
+70.62 0 l
+0 16.29 l
+47.46 0 l
+0 24.43 l
+38.02 0 l
+0 28.51 l
+380.52 0 l
+0 85.53 l
+0 0 l
+o
+0.75 setlinewidth
+[ 3.00 5.00] 0 setdash
+np
+103.58 107.25 m
+12.37 0 l
+0 11.62 l
+23.89 0 l
+0 12.40 l
+6.12 0 l
+0 12.40 l
+10.71 0 l
+0 12.39 l
+30.49 0 l
+0 12.40 l
+2.48 0 l
+0 12.40 l
+19.12 0 l
+0 12.40 l
+87.89 0 l
+0 17.35 l
+70.00 0 l
+0 19.84 l
+115.86 0 l
+0 19.84 l
+0 0 l
+o
+0.75 setlinewidth
+[ 0.00 4.00] 0 setdash
+np
+103.58 107.25 m
+7.90 0 l
+0 13.62 l
+15.35 0 l
+0 13.63 l
+10.22 0 l
+0 14.22 l
+16.77 0 l
+0 14.22 l
+3.61 0 l
+0 14.22 l
+11.28 0 l
+0 15.06 l
+19.10 0 l
+0 16.05 l
+2.53 0 l
+0 16.06 l
+4.74 0 l
+0 16.06 l
+34.27 0 l
+0 17.67 l
+44.52 0 l
+0 17.66 l
+64.18 0 l
+0 24.74 l
+46.88 0 l
+0 32.97 l
+145.45 0 l
+o
+0.75 setlinewidth
+[ 0.00 4.00 3.00 4.00] 0 setdash
+np
+103.58 107.25 m
+14.84 0 l
+0 11.97 l
+4.04 0 l
+0 11.98 l
+0.06 0 l
+0 11.97 l
+41.71 0 l
+0 12.87 l
+1.77 0 l
+0 12.86 l
+28.70 0 l
+0 13.44 l
+17.93 0 l
+0 14.12 l
+27.14 0 l
+0 14.12 l
+0.34 0 l
+0 14.12 l
+169.63 0 l
+0 28.51 l
+0.35 0 l
+0 28.51 l
+357.02 0 l
+o
+0.75 setlinewidth
+[ 6.00 4.00] 0 setdash
+np
+103.58 107.25 m
+9.19 0 l
+0 11.62 l
+29.82 0 l
+0 12.40 l
+15.66 0 l
+0 12.40 l
+4.36 0 l
+0 12.39 l
+8.76 0 l
+0 12.40 l
+20.40 0 l
+0 12.40 l
+15.07 0 l
+0 12.40 l
+7.65 0 l
+0 12.40 l
+5.14 0 l
+0 12.39 l
+6.68 0 l
+0 12.40 l
+6.52 0 l
+0 12.40 l
+67.08 0 l
+0 17.36 l
+12.72 0 l
+0 19.83 l
+6.50 0 l
+0 19.84 l
+10.27 0 l
+0 19.84 l
+17.09 0 l
+0 19.83 l
+77.59 0 l
+0 19.84 l
+58.43 0 l
+o
+0.75 setlinewidth
+[ 1.00 3.00 5.00 3.00] 0 setdash
+np
+103.58 107.25 m
+5.46 0 l
+0 13.62 l
+8.95 0 l
+0 13.63 l
+7.65 0 l
+0 13.63 l
+10.26 0 l
+0 14.22 l
+0.47 0 l
+0 14.22 l
+64.95 0 l
+0 16.06 l
+35.63 0 l
+0 17.66 l
+284.30 0 l
+0 32.98 l
+9.13 0 l
+o
+ep
+%%Page: 2 2
+bp
+77.04 91.44 793.65 518.24 cl
+18.00 18.00 823.89 577.28 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 91.44 m
+632.09 0 l
+o
+np
+103.58 91.44 m
+0 -7.20 l
+o
+np
+230.00 91.44 m
+0 -7.20 l
+o
+np
+356.42 91.44 m
+0 -7.20 l
+o
+np
+482.84 91.44 m
+0 -7.20 l
+o
+np
+609.26 91.44 m
+0 -7.20 l
+o
+np
+735.67 91.44 m
+0 -7.20 l
+o
+/ps 12 def /Font1 findfont 12 s
+103.58 65.52 (0) .5 0 0 t
+230.00 65.52 (1) .5 0 0 t
+356.42 65.52 (2) .5 0 0 t
+482.84 65.52 (3) .5 0 0 t
+609.26 65.52 (4) .5 0 0 t
+735.67 65.52 (5) .5 0 0 t
+np
+77.04 107.25 m
+0 395.18 l
+o
+np
+77.04 107.25 m
+-7.20 0 l
+o
+np
+77.04 186.28 m
+-7.20 0 l
+o
+np
+77.04 265.32 m
+-7.20 0 l
+o
+np
+77.04 344.36 m
+-7.20 0 l
+o
+np
+77.04 423.39 m
+-7.20 0 l
+o
+np
+77.04 502.43 m
+-7.20 0 l
+o
+59.76 107.25 (0.0) .5 0 90 t
+59.76 186.28 (0.2) .5 0 90 t
+59.76 265.32 (0.4) .5 0 90 t
+59.76 344.36 (0.6) .5 0 90 t
+59.76 423.39 (0.8) .5 0 90 t
+59.76 502.43 (1.0) .5 0 90 t
+np
+77.04 518.24 m
+0 -426.80 l
+716.61 0 l
+o
+18.00 18.00 823.89 577.28 cl
+/ps 14 def /Font2 findfont 14 s
+0 setgray
+435.34 542.73 ( ) .5 0 0 t
+/ps 12 def /Font1 findfont 12 s
+435.34 36.72 (Years) .5 0 0 t
+30.96 304.84 (Probability) .5 0 90 t
+77.04 91.44 793.65 518.24 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+114.38 488.03 m
+21.60 0 l
+o
+1 0 0 rgb
+np
+114.38 473.63 m
+21.60 0 l
+o
+0 0.8039 0 rgb
+np
+114.38 459.23 m
+21.60 0 l
+o
+0 0 1 rgb
+np
+114.38 444.83 m
+21.60 0 l
+o
+0 1 1 rgb
+np
+114.38 430.43 m
+21.60 0 l
+o
+1 0 1 rgb
+np
+114.38 416.03 m
+21.60 0 l
+o
+/ps 12 def /Font1 findfont 12 s
+0 setgray
+146.78 483.72 (a 1) 0 0 0 t
+146.78 469.32 (b 1) 0 0 0 t
+146.78 454.92 (c 1) 0 0 0 t
+146.78 440.52 (a 2) 0 0 0 t
+146.78 426.12 (b 2) 0 0 0 t
+146.78 411.72 (c 2) 0 0 0 t
+np
+103.58 107.25 m
+17.20 0 l
+0 11.97 l
+22.03 0 l
+0 12.86 l
+26.56 0 l
+0 12.87 l
+30.66 0 l
+0 14.12 l
+30.46 0 l
+0 14.12 l
+70.62 0 l
+0 16.29 l
+47.46 0 l
+0 24.43 l
+38.02 0 l
+0 28.51 l
+380.52 0 l
+0 85.53 l
+0 0 l
+o
+1 0 0 rgb
+np
+103.58 107.25 m
+12.37 0 l
+0 11.62 l
+23.89 0 l
+0 12.40 l
+6.12 0 l
+0 12.40 l
+10.71 0 l
+0 12.39 l
+30.49 0 l
+0 12.40 l
+2.48 0 l
+0 12.40 l
+19.12 0 l
+0 12.40 l
+87.89 0 l
+0 17.35 l
+70.00 0 l
+0 19.84 l
+115.86 0 l
+0 19.84 l
+0 0 l
+o
+0 0.8039 0 rgb
+np
+103.58 107.25 m
+7.90 0 l
+0 13.62 l
+15.35 0 l
+0 13.63 l
+10.22 0 l
+0 14.22 l
+16.77 0 l
+0 14.22 l
+3.61 0 l
+0 14.22 l
+11.28 0 l
+0 15.06 l
+19.10 0 l
+0 16.05 l
+2.53 0 l
+0 16.06 l
+4.74 0 l
+0 16.06 l
+34.27 0 l
+0 17.67 l
+44.52 0 l
+0 17.66 l
+64.18 0 l
+0 24.74 l
+46.88 0 l
+0 32.97 l
+145.45 0 l
+o
+0 0 1 rgb
+np
+103.58 107.25 m
+14.84 0 l
+0 11.97 l
+4.04 0 l
+0 11.98 l
+0.06 0 l
+0 11.97 l
+41.71 0 l
+0 12.87 l
+1.77 0 l
+0 12.86 l
+28.70 0 l
+0 13.44 l
+17.93 0 l
+0 14.12 l
+27.14 0 l
+0 14.12 l
+0.34 0 l
+0 14.12 l
+169.63 0 l
+0 28.51 l
+0.35 0 l
+0 28.51 l
+357.02 0 l
+o
+0 1 1 rgb
+np
+103.58 107.25 m
+9.19 0 l
+0 11.62 l
+29.82 0 l
+0 12.40 l
+15.66 0 l
+0 12.40 l
+4.36 0 l
+0 12.39 l
+8.76 0 l
+0 12.40 l
+20.40 0 l
+0 12.40 l
+15.07 0 l
+0 12.40 l
+7.65 0 l
+0 12.40 l
+5.14 0 l
+0 12.39 l
+6.68 0 l
+0 12.40 l
+6.52 0 l
+0 12.40 l
+67.08 0 l
+0 17.36 l
+12.72 0 l
+0 19.83 l
+6.50 0 l
+0 19.84 l
+10.27 0 l
+0 19.84 l
+17.09 0 l
+0 19.83 l
+77.59 0 l
+0 19.84 l
+58.43 0 l
+o
+1 0 1 rgb
+np
+103.58 107.25 m
+5.46 0 l
+0 13.62 l
+8.95 0 l
+0 13.63 l
+7.65 0 l
+0 13.63 l
+10.26 0 l
+0 14.22 l
+0.47 0 l
+0 14.22 l
+64.95 0 l
+0 16.06 l
+35.63 0 l
+0 17.66 l
+284.30 0 l
+0 32.98 l
+9.13 0 l
+o
+ep
+%%Page: 3 3
+bp
+77.04 91.44 793.65 518.24 cl
+77.04 91.44 793.65 518.24 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+103.58 107.25 2.70 c p1
+108.11 109.52 2.70 c p1
+112.99 109.97 2.70 c p1
+119.12 451.09 2.70 c p1
+129.46 463.64 2.70 c p1
+132.28 467.06 2.70 c p1
+135.29 469.74 2.70 c p1
+138.47 131.03 2.70 c p1
+146.44 473.16 2.70 c p1
+149.32 476.64 2.70 c p1
+150.09 137.12 2.70 c p1
+161.50 483.65 2.70 c p1
+162.17 143.23 2.70 c p1
+180.18 488.29 2.70 c p1
+180.84 148.00 2.70 c p1
+182.68 487.08 2.70 c p1
+188.19 149.16 2.70 c p1
+193.20 493.02 2.70 c p1
+202.03 496.99 2.70 c p1
+222.87 162.78 2.70 c p1
+224.02 502.43 2.70 c p1
+267.93 160.40 2.70 c p1
+290.99 155.10 2.70 c p1
+295.50 493.45 2.70 c p1
+332.88 148.39 2.70 c p1
+343.53 132.35 2.70 c p1
+361.83 451.61 2.70 c p1
+380.32 117.41 2.70 c p1
+382.01 450.84 2.70 c p1
+479.08 412.19 2.70 c p1
+767.11 317.57 2.70 c p1
+18.00 18.00 823.89 577.28 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+95.59 91.44 m
+639.71 0 l
+o
+np
+95.59 91.44 m
+0 -7.20 l
+o
+np
+223.53 91.44 m
+0 -7.20 l
+o
+np
+351.47 91.44 m
+0 -7.20 l
+o
+np
+479.41 91.44 m
+0 -7.20 l
+o
+np
+607.35 91.44 m
+0 -7.20 l
+o
+np
+735.30 91.44 m
+0 -7.20 l
+o
+/ps 12 def /Font1 findfont 12 s
+95.59 65.52 (0) .5 0 0 t
+223.53 65.52 (1) .5 0 0 t
+351.47 65.52 (2) .5 0 0 t
+479.41 65.52 (3) .5 0 0 t
+607.35 65.52 (4) .5 0 0 t
+735.30 65.52 (5) .5 0 0 t
+np
+77.04 110.48 m
+0 342.22 l
+o
+np
+77.04 110.48 m
+-7.20 0 l
+o
+np
+77.04 178.93 m
+-7.20 0 l
+o
+np
+77.04 247.37 m
+-7.20 0 l
+o
+np
+77.04 315.81 m
+-7.20 0 l
+o
+np
+77.04 384.25 m
+-7.20 0 l
+o
+np
+77.04 452.70 m
+-7.20 0 l
+o
+59.76 110.48 (-0.6) .5 0 90 t
+59.76 178.93 (-0.4) .5 0 90 t
+59.76 247.37 (-0.2) .5 0 90 t
+59.76 315.81 (0.0) .5 0 90 t
+59.76 384.25 (0.2) .5 0 90 t
+59.76 452.70 (0.4) .5 0 90 t
+np
+77.04 91.44 m
+716.61 0 l
+0 426.80 l
+-716.61 0 l
+0 -426.80 l
+o
+18.00 18.00 823.89 577.28 cl
+/ps 12 def /Font1 findfont 12 s
+0 setgray
+435.34 36.72 (ww$uft) .5 0 0 t
+30.96 304.84 (ww$res[, 1]) .5 0 90 t
+77.04 91.44 793.65 518.24 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 258.61 m
+4.53 6.89 l
+4.88 7.24 l
+6.13 8.83 l
+10.34 14.49 l
+2.82 3.88 l
+3.01 4.13 l
+3.18 4.34 l
+7.97 10.66 l
+2.88 3.76 l
+0.77 1.01 l
+11.41 14.03 l
+0.67 0.76 l
+18.01 15.09 l
+0.66 0.21 l
+1.84 0.61 l
+5.51 0.21 l
+5.01 -1.46 l
+8.83 -5.73 l
+20.84 -22.08 l
+1.15 -0.50 l
+43.91 -32.84 l
+23.06 -2.72 l
+4.51 -1.08 l
+37.38 -4.40 l
+10.65 3.35 l
+18.30 6.57 l
+18.49 6.83 l
+1.69 0.63 l
+97.07 44.04 l
+288.03 -3.56 l
+o
+ep
+%%Page: 4 4
+bp
+77.04 91.44 793.65 518.24 cl
+18.00 18.00 823.89 577.28 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 91.44 m
+632.09 0 l
+o
+np
+103.58 91.44 m
+0 -7.20 l
+o
+np
+230.00 91.44 m
+0 -7.20 l
+o
+np
+356.42 91.44 m
+0 -7.20 l
+o
+np
+482.84 91.44 m
+0 -7.20 l
+o
+np
+609.26 91.44 m
+0 -7.20 l
+o
+np
+735.67 91.44 m
+0 -7.20 l
+o
+/ps 12 def /Font1 findfont 12 s
+103.58 65.52 (0) .5 0 0 t
+230.00 65.52 (1) .5 0 0 t
+356.42 65.52 (2) .5 0 0 t
+482.84 65.52 (3) .5 0 0 t
+609.26 65.52 (4) .5 0 0 t
+735.67 65.52 (5) .5 0 0 t
+np
+77.04 107.25 m
+0 410.36 l
+o
+np
+77.04 107.25 m
+-7.20 0 l
+o
+np
+77.04 175.64 m
+-7.20 0 l
+o
+np
+77.04 244.04 m
+-7.20 0 l
+o
+np
+77.04 312.43 m
+-7.20 0 l
+o
+np
+77.04 380.83 m
+-7.20 0 l
+o
+np
+77.04 449.22 m
+-7.20 0 l
+o
+np
+77.04 517.61 m
+-7.20 0 l
+o
+59.76 107.25 (0.0) .5 0 90 t
+59.76 175.64 (0.1) .5 0 90 t
+59.76 244.04 (0.2) .5 0 90 t
+59.76 312.43 (0.3) .5 0 90 t
+59.76 380.83 (0.4) .5 0 90 t
+59.76 449.22 (0.5) .5 0 90 t
+59.76 517.61 (0.6) .5 0 90 t
+np
+77.04 91.44 m
+716.61 0 l
+0 426.80 l
+-716.61 0 l
+0 -426.80 l
+o
+18.00 18.00 823.89 577.28 cl
+/ps 12 def /Font1 findfont 12 s
+0 setgray
+435.34 36.72 (c\(xmin, xmax\)) .5 0 0 t
+30.96 304.84 (ylim) .5 0 90 t
+77.04 91.44 793.65 518.24 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 107.25 m
+7.90 0 l
+0 5.70 l
+4.47 0 l
+0 5.73 l
+4.83 0 l
+0 5.67 l
+6.05 0 l
+0 5.65 l
+10.22 0 l
+0 5.80 l
+2.79 0 l
+0 5.80 l
+2.97 0 l
+0 5.75 l
+3.15 0 l
+0 5.74 l
+7.86 0 l
+0 5.69 l
+2.85 0 l
+0 5.68 l
+0.76 0 l
+0 5.71 l
+11.28 0 l
+0 5.69 l
+0.66 0 l
+0 5.68 l
+17.79 0 l
+0 5.68 l
+0.65 0 l
+0 5.67 l
+1.83 0 l
+0 5.72 l
+5.44 0 l
+0 5.71 l
+4.95 0 l
+0 5.80 l
+8.73 0 l
+0 5.73 l
+20.59 0 l
+0 5.67 l
+1.14 0 l
+0 5.71 l
+43.38 0 l
+0 5.88 l
+22.78 0 l
+0 6.96 l
+4.46 0 l
+0 7.04 l
+36.94 0 l
+0 8.99 l
+10.52 0 l
+0 9.89 l
+18.08 0 l
+0 11.18 l
+18.28 0 l
+0 13.06 l
+1.66 0 l
+0 13.29 l
+95.92 0 l
+0 19.10 l
+284.60 0 l
+0 50.71 l
+0 0 l
+o
+0.75 setlinewidth
+[ 3.00 5.00] 0 setdash
+np
+103.58 107.25 m
+7.90 0 l
+0 9.44 l
+4.47 0 l
+0 9.70 l
+4.83 0 l
+0 9.82 l
+6.05 0 l
+0 10.09 l
+10.22 0 l
+0 10.88 l
+2.79 0 l
+0 10.95 l
+2.97 0 l
+0 10.94 l
+3.15 0 l
+0 10.99 l
+7.86 0 l
+0 11.23 l
+2.85 0 l
+0 11.26 l
+0.76 0 l
+0 11.23 l
+11.28 0 l
+0 11.65 l
+0.66 0 l
+0 11.52 l
+17.79 0 l
+0 12.20 l
+0.65 0 l
+0 12.05 l
+1.83 0 l
+0 12.07 l
+5.44 0 l
+0 12.09 l
+4.95 0 l
+0 12.31 l
+8.73 0 l
+0 12.26 l
+20.59 0 l
+0 12.56 l
+1.14 0 l
+0 12.45 l
+43.38 0 l
+0 13.24 l
+22.78 0 l
+0 15.31 l
+4.46 0 l
+0 15.08 l
+36.94 0 l
+0 17.83 l
+10.52 0 l
+0 18.61 l
+18.08 0 l
+0 19.37 l
+18.28 0 l
+0 20.52 l
+1.66 0 l
+0 20.02 l
+95.92 0 l
+0 16.31 l
+284.60 0 l
+0 1.20 l
+0 0 l
+o
+ep
+%%Page: 5 5
+bp
+77.04 91.44 793.65 518.24 cl
+18.00 18.00 823.89 577.28 cl
+0 setgray
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 91.44 m
+632.09 0 l
+o
+np
+103.58 91.44 m
+0 -7.20 l
+o
+np
+230.00 91.44 m
+0 -7.20 l
+o
+np
+356.42 91.44 m
+0 -7.20 l
+o
+np
+482.84 91.44 m
+0 -7.20 l
+o
+np
+609.26 91.44 m
+0 -7.20 l
+o
+np
+735.67 91.44 m
+0 -7.20 l
+o
+/ps 12 def /Font1 findfont 12 s
+103.58 65.52 (0) .5 0 0 t
+230.00 65.52 (1) .5 0 0 t
+356.42 65.52 (2) .5 0 0 t
+482.84 65.52 (3) .5 0 0 t
+609.26 65.52 (4) .5 0 0 t
+735.67 65.52 (5) .5 0 0 t
+np
+77.04 107.25 m
+0 410.36 l
+o
+np
+77.04 107.25 m
+-7.20 0 l
+o
+np
+77.04 175.64 m
+-7.20 0 l
+o
+np
+77.04 244.04 m
+-7.20 0 l
+o
+np
+77.04 312.43 m
+-7.20 0 l
+o
+np
+77.04 380.83 m
+-7.20 0 l
+o
+np
+77.04 449.22 m
+-7.20 0 l
+o
+np
+77.04 517.61 m
+-7.20 0 l
+o
+59.76 107.25 (0.0) .5 0 90 t
+59.76 175.64 (0.1) .5 0 90 t
+59.76 244.04 (0.2) .5 0 90 t
+59.76 312.43 (0.3) .5 0 90 t
+59.76 380.83 (0.4) .5 0 90 t
+59.76 449.22 (0.5) .5 0 90 t
+59.76 517.61 (0.6) .5 0 90 t
+np
+77.04 91.44 m
+716.61 0 l
+0 426.80 l
+-716.61 0 l
+0 -426.80 l
+o
+18.00 18.00 823.89 577.28 cl
+/ps 12 def /Font1 findfont 12 s
+0 setgray
+435.34 36.72 (c\(xmin, xmax\)) .5 0 0 t
+30.96 304.84 (ylim) .5 0 90 t
+77.04 91.44 793.65 518.24 cl
+1 0 0 rgb
+0.75 setlinewidth
+[] 0 setdash
+1 setlinecap
+1 setlinejoin
+10.00 setmiterlimit
+np
+103.58 107.25 m
+7.90 0 l
+0 5.70 l
+4.47 0 l
+0 5.73 l
+4.83 0 l
+0 5.67 l
+6.05 0 l
+0 5.65 l
+10.22 0 l
+0 5.80 l
+2.79 0 l
+0 5.80 l
+2.97 0 l
+0 5.75 l
+3.15 0 l
+0 5.74 l
+7.86 0 l
+0 5.69 l
+2.85 0 l
+0 5.68 l
+0.76 0 l
+0 5.71 l
+11.28 0 l
+0 5.69 l
+0.66 0 l
+0 5.68 l
+17.79 0 l
+0 5.68 l
+0.65 0 l
+0 5.67 l
+1.83 0 l
+0 5.72 l
+5.44 0 l
+0 5.71 l
+4.95 0 l
+0 5.80 l
+8.73 0 l
+0 5.73 l
+20.59 0 l
+0 5.67 l
+1.14 0 l
+0 5.71 l
+43.38 0 l
+0 5.88 l
+22.78 0 l
+0 6.96 l
+4.46 0 l
+0 7.04 l
+36.94 0 l
+0 8.99 l
+10.52 0 l
+0 9.89 l
+18.08 0 l
+0 11.18 l
+18.28 0 l
+0 13.06 l
+1.66 0 l
+0 13.29 l
+95.92 0 l
+0 19.10 l
+284.60 0 l
+0 50.71 l
+0 0 l
+o
+0 0 1 rgb
+np
+103.58 107.25 m
+7.90 0 l
+0 9.44 l
+4.47 0 l
+0 9.70 l
+4.83 0 l
+0 9.82 l
+6.05 0 l
+0 10.09 l
+10.22 0 l
+0 10.88 l
+2.79 0 l
+0 10.95 l
+2.97 0 l
+0 10.94 l
+3.15 0 l
+0 10.99 l
+7.86 0 l
+0 11.23 l
+2.85 0 l
+0 11.26 l
+0.76 0 l
+0 11.23 l
+11.28 0 l
+0 11.65 l
+0.66 0 l
+0 11.52 l
+17.79 0 l
+0 12.20 l
+0.65 0 l
+0 12.05 l
+1.83 0 l
+0 12.07 l
+5.44 0 l
+0 12.09 l
+4.95 0 l
+0 12.31 l
+8.73 0 l
+0 12.26 l
+20.59 0 l
+0 12.56 l
+1.14 0 l
+0 12.45 l
+43.38 0 l
+0 13.24 l
+22.78 0 l
+0 15.31 l
+4.46 0 l
+0 15.08 l
+36.94 0 l
+0 17.83 l
+10.52 0 l
+0 18.61 l
+18.08 0 l
+0 19.37 l
+18.28 0 l
+0 20.52 l
+1.66 0 l
+0 20.02 l
+95.92 0 l
+0 16.31 l
+284.60 0 l
+0 1.20 l
+0 0 l
+o
+ep
+%%Trailer
+%%Pages: 5
+%%EOF
diff --git a/tests/test.R b/tests/test.R
new file mode 100644
index 0000000..f9a5f1c
--- /dev/null
+++ b/tests/test.R
@@ -0,0 +1,53 @@
+library(cmprsk)
+u <- R.Version()
+if (u$major>'1' | u$major == '1' & u$minor>'6.2') RNGversion("1.6.2")
+
+set.seed(2)
+ss <- rexp(100)
+gg <- factor(sample(1:3,100,replace=TRUE),1:3,c('a','b','c'))
+cc <- sample(0:2,100,replace=TRUE)
+strt <- sample(1:2,100,replace=TRUE)
+dd <- data.frame(ss=abs(rnorm(100)))
+d2 <- data.frame(ssd=ss,ggd=gg,ccd=cc,strtd=strt,X=c(rep(1,80),rep(0,20)))
+gg2 <- gg
+gg2[c(5,10,50)] <- NA
+print(xx <- cuminc(dd$ss,cc,gg,strt))
+print(xx <- cuminc(d2$ssd,d2$ccd))
+print(xx <- cuminc(ss,cc))
+print(xx <- cuminc(ss,cc,gg,strt))
+plot(xx)
+plot(xx,lty=1,color=1:6)
+print(xx <- cuminc(dd$ss,cc,gg2,strt))
+print(xx <- cuminc(ss,cc,gg2,strt,subset=d2$X == 1))
+attach(d2)
+print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=X == 1))
+print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=gg != 'b'))
+print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=ggd != 'b'))
+print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=gg2 != 'b'))
+detach(d2)
+
+cv <- matrix(sample(0:1,3*100,replace=TRUE),ncol=3)
+cv[c(1,10,20)] <- NA
+print(xx <- crr(ss,cc,cv))
+cov2 <- cbind(cv[,1],cv[,1])
+tf <- function(uft) cbind(uft,uft^2)
+print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3]))
+plot(ww$uft,ww$res[,1])
+lines(lowess(ww$uft,ww$res[,1],iter=0,f=.75))
+print(wp <- predict(ww,rbind(c(1,1,1),c(0,0,0)),rbind(c(1,1),c(0,0))))
+plot(wp)
+plot(wp,lty=1,col=c(2,4))
+d2 <- cbind(d2,cv3=cv[,3],cv1=cv[,1])
+attach(d2)
+print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cov2,tf=tf,cengroup=cv3))
+print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cov2,tf=tf))
+print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cbind(cv1,cv1),tf=tf,cengroup=cv3))
+print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cbind(cv1,cv1),tf=tf,cengroup=cv3,subset=X == 1))
+print(summary(ww))
+detach(d2)
+print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3],subset=d2$X==1))
+print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3],failcode=2))
+print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3],cencode=2))
+print(ww <- crr(ss,cc,cv[,1]))
+print(ww <- crr(ss,cc,cov2=cv[,1],tf=function(x) x))
+print(summary(ww))
diff --git a/tests/test.Rout.save b/tests/test.Rout.save
new file mode 100644
index 0000000..9aadd6b
--- /dev/null
+++ b/tests/test.Rout.save
@@ -0,0 +1,476 @@
+
+R version 2.15.0 (2012-03-30)
+Copyright (C) 2012 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(cmprsk)
+Loading required package: survival
+Loading required package: splines
+> u <- R.Version()
+> if (u$major>'1' | u$major == '1' & u$minor>'6.2') RNGversion("1.6.2")
+Warning message:
+In RNGkind("Marsaglia-Multicarry", "Buggy Kinderman-Ramage") :
+ Buggy version of Kinderman-Ramage generator used
+>
+> set.seed(2)
+> ss <- rexp(100)
+> gg <- factor(sample(1:3,100,replace=TRUE),1:3,c('a','b','c'))
+> cc <- sample(0:2,100,replace=TRUE)
+> strt <- sample(1:2,100,replace=TRUE)
+> dd <- data.frame(ss=abs(rnorm(100)))
+> d2 <- data.frame(ssd=ss,ggd=gg,ccd=cc,strtd=strt,X=c(rep(1,80),rep(0,20)))
+> gg2 <- gg
+> gg2[c(5,10,50)] <- NA
+> print(xx <- cuminc(dd$ss,cc,gg,strt))
+Tests:
+ stat pv df
+1 1.214132 0.5449473 2
+2 3.269462 0.1950048 2
+Estimates and Variances:
+$est
+ 1 2 3
+a 1 0.2105985 NA NA
+b 1 0.2685413 0.3471880 0.3471880
+c 1 0.3082577 0.4989848 NA
+a 2 0.3156040 NA NA
+b 2 0.3382252 0.6036578 0.6036578
+c 2 0.2376301 0.2830413 NA
+
+$var
+ 1 2 3
+a 1 0.006386094 NA NA
+b 1 0.007003493 0.008623542 0.008623542
+c 1 0.008874672 0.012216283 NA
+a 2 0.008549829 NA NA
+b 2 0.008107166 0.010928915 0.010928915
+c 2 0.007656447 0.008944878 NA
+
+> print(xx <- cuminc(d2$ssd,d2$ccd))
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+1 1 0.2352861 0.3287068 0.4247898 0.4247898 0.4247898
+1 2 0.2567862 0.3931580 0.4631781 0.4911861 0.4911861
+
+$var
+ 1 2 3 4 5
+1 1 0.002062112 0.002957844 0.004170634 0.004170634 0.004170634
+1 2 0.002190305 0.003273248 0.004004514 0.004369940 0.004369940
+
+> print(xx <- cuminc(ss,cc))
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+1 1 0.2352861 0.3287068 0.4247898 0.4247898 0.4247898
+1 2 0.2567862 0.3931580 0.4631781 0.4911861 0.4911861
+
+$var
+ 1 2 3 4 5
+1 1 0.002062112 0.002957844 0.004170634 0.004170634 0.004170634
+1 2 0.002190305 0.003273248 0.004004514 0.004369940 0.004369940
+
+> print(xx <- cuminc(ss,cc,gg,strt))
+Tests:
+ stat pv df
+1 3.393977 0.1832345 2
+2 1.989511 0.3698139 2
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+a 1 0.1311269 0.2699184 0.3420625 0.3420625 0.3420625
+b 1 0.2176471 0.2615686 NA NA NA
+c 1 0.3816280 0.4889137 0.5723581 NA NA
+a 2 0.2257601 0.2972171 0.4415053 0.4415053 0.4415053
+b 2 0.3117647 0.5878431 NA NA NA
+c 2 0.2160508 0.2607532 0.2607532 NA NA
+
+$var
+ 1 2 3 4 5
+a 1 0.003922836 0.009113464 0.012507959 0.01250796 0.01250796
+b 1 0.005528913 0.006916070 NA NA NA
+c 1 0.009854126 0.012351948 0.015575726 NA NA
+a 2 0.005947601 0.007356628 0.014479506 0.01447951 0.01447951
+b 2 0.007033454 0.010927780 NA NA NA
+c 2 0.006497566 0.007867543 0.007867543 NA NA
+
+> plot(xx)
+> plot(xx,lty=1,color=1:6)
+> print(xx <- cuminc(dd$ss,cc,gg2,strt))
+3 cases omitted due to missing values
+Tests:
+ stat pv df
+1 1.008453 0.6039725 2
+2 2.728735 0.2555423 2
+Estimates and Variances:
+$est
+ 1 2 3
+a 1 0.2105985 NA NA
+b 1 0.2987226 0.3888341 0.3888341
+c 1 0.3082577 0.4989848 NA
+a 2 0.3156040 NA NA
+b 2 0.3408314 0.5661101 0.5661101
+c 2 0.2376301 0.2830413 NA
+
+$var
+ 1 2 3
+a 1 0.006386094 NA NA
+b 1 0.008392163 0.010391194 0.01039119
+c 1 0.008874672 0.012216283 NA
+a 2 0.008549829 NA NA
+b 2 0.009203557 0.011899088 0.01189909
+c 2 0.007656447 0.008944878 NA
+
+> print(xx <- cuminc(ss,cc,gg2,strt,subset=d2$X == 1))
+3 cases omitted due to missing values
+Tests:
+ stat pv df
+1 4.716769 0.0945729 2
+2 2.696396 0.2597078 2
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+a 1 0.1279762 0.1848443 0.1848443 0.1848443 0.1848443
+b 1 0.1431818 0.2147727 NA NA NA
+c 1 0.3740171 0.4874823 0.5757330 NA NA
+a 2 0.2085623 0.3033425 0.5080678 0.5080678 0.5080678
+b 2 0.3795455 0.6062500 NA NA NA
+c 2 0.2004884 0.2477656 0.2477656 NA NA
+
+$var
+ 1 2 3 4 5
+a 1 0.005087085 0.007761560 0.007761560 0.00776156 0.00776156
+b 1 0.006224509 0.010655685 NA NA NA
+c 1 0.010614083 0.013441268 0.017054461 NA NA
+a 2 0.007382690 0.009967497 0.022418302 0.02241830 0.02241830
+b 2 0.012049806 0.019716303 NA NA NA
+c 2 0.006846162 0.008398669 0.008398669 NA NA
+
+> attach(d2)
+> print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=X == 1))
+3 cases omitted due to missing values
+Tests:
+ stat pv df
+1 4.716769 0.0945729 2
+2 2.696396 0.2597078 2
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+a 1 0.1279762 0.1848443 0.1848443 0.1848443 0.1848443
+b 1 0.1431818 0.2147727 NA NA NA
+c 1 0.3740171 0.4874823 0.5757330 NA NA
+a 2 0.2085623 0.3033425 0.5080678 0.5080678 0.5080678
+b 2 0.3795455 0.6062500 NA NA NA
+c 2 0.2004884 0.2477656 0.2477656 NA NA
+
+$var
+ 1 2 3 4 5
+a 1 0.005087085 0.007761560 0.007761560 0.00776156 0.00776156
+b 1 0.006224509 0.010655685 NA NA NA
+c 1 0.010614083 0.013441268 0.017054461 NA NA
+a 2 0.007382690 0.009967497 0.022418302 0.02241830 0.02241830
+b 2 0.012049806 0.019716303 NA NA NA
+c 2 0.006846162 0.008398669 0.008398669 NA NA
+
+> print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=gg != 'b'))
+Tests:
+ stat pv df
+1 3.5151549 0.06080997 1
+2 0.1400145 0.70826655 1
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+a 1 0.1311269 0.2699184 0.3420625 0.3420625 0.3420625
+c 1 0.3816280 0.4889137 0.5723581 NA NA
+a 2 0.2257601 0.2972171 0.4415053 0.4415053 0.4415053
+c 2 0.2160508 0.2607532 0.2607532 NA NA
+
+$var
+ 1 2 3 4 5
+a 1 0.003922836 0.009113464 0.012507959 0.01250796 0.01250796
+c 1 0.009854126 0.012351948 0.015575726 NA NA
+a 2 0.005947601 0.007356628 0.014479506 0.01447951 0.01447951
+c 2 0.006497566 0.007867543 0.007867543 NA NA
+
+> print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=ggd != 'b'))
+Tests:
+ stat pv df
+1 3.5151549 0.06080997 1
+2 0.1400145 0.70826655 1
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+a 1 0.1311269 0.2699184 0.3420625 0.3420625 0.3420625
+c 1 0.3816280 0.4889137 0.5723581 NA NA
+a 2 0.2257601 0.2972171 0.4415053 0.4415053 0.4415053
+c 2 0.2160508 0.2607532 0.2607532 NA NA
+
+$var
+ 1 2 3 4 5
+a 1 0.003922836 0.009113464 0.012507959 0.01250796 0.01250796
+c 1 0.009854126 0.012351948 0.015575726 NA NA
+a 2 0.005947601 0.007356628 0.014479506 0.01447951 0.01447951
+c 2 0.006497566 0.007867543 0.007867543 NA NA
+
+> print(xx <- cuminc(ssd,ccd,gg2,strtd,subset=gg2 != 'b'))
+3 cases omitted due to missing values
+Tests:
+ stat pv df
+1 3.5151549 0.06080997 1
+2 0.1400145 0.70826655 1
+Estimates and Variances:
+$est
+ 1 2 3 4 5
+a 1 0.1311269 0.2699184 0.3420625 0.3420625 0.3420625
+c 1 0.3816280 0.4889137 0.5723581 NA NA
+a 2 0.2257601 0.2972171 0.4415053 0.4415053 0.4415053
+c 2 0.2160508 0.2607532 0.2607532 NA NA
+
+$var
+ 1 2 3 4 5
+a 1 0.003922836 0.009113464 0.012507959 0.01250796 0.01250796
+c 1 0.009854126 0.012351948 0.015575726 NA NA
+a 2 0.005947601 0.007356628 0.014479506 0.01447951 0.01447951
+c 2 0.006497566 0.007867543 0.007867543 NA NA
+
+> detach(d2)
+>
+> cv <- matrix(sample(0:1,3*100,replace=TRUE),ncol=3)
+> cv[c(1,10,20)] <- NA
+> print(xx <- crr(ss,cc,cv))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+ cv1 cv2 cv3
+-0.29130 -0.06346 -0.47510
+standard errors:
+[1] 0.3698 0.3566 0.3544
+two-sided p-values:
+ cv1 cv2 cv3
+0.43 0.86 0.18
+> cov2 <- cbind(cv[,1],cv[,1])
+> tf <- function(uft) cbind(uft,uft^2)
+> print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3]))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.02944 -0.06263 -0.42160 -0.84780 0.29850
+standard errors:
+[1] 0.8029 0.3583 0.3595 1.3090 0.4004
+two-sided p-values:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.97 0.86 0.24 0.52 0.46
+> plot(ww$uft,ww$res[,1])
+> lines(lowess(ww$uft,ww$res[,1],iter=0,f=.75))
+> print(wp <- predict(ww,rbind(c(1,1,1),c(0,0,0)),rbind(c(1,1),c(0,0))))
+ [,1] [,2] [,3]
+ [1,] 0.06246194 0.008345203 0.01381136
+ [2,] 0.09786020 0.016721921 0.02799111
+ [3,] 0.13602406 0.025009866 0.04235389
+ [4,] 0.18390522 0.033273945 0.05710125
+ [5,] 0.26471120 0.041754133 0.07301143
+ [6,] 0.28678273 0.050225177 0.08901147
+ [7,] 0.31027641 0.058637012 0.10501213
+ [8,] 0.33519487 0.067025773 0.12108626
+ [9,] 0.39744032 0.075343709 0.13750391
+[10,] 0.41997041 0.083655640 0.15396776
+[11,] 0.42598719 0.092000584 0.17038559
+[12,] 0.51517185 0.100319016 0.18741181
+[13,] 0.52041769 0.108625866 0.20426566
+[14,] 0.66116579 0.116924371 0.22209367
+[15,] 0.66629928 0.125218167 0.23971535
+[16,] 0.68071685 0.133587124 0.25736827
+[17,] 0.72376095 0.141928165 0.27503823
+[18,] 0.76292157 0.150413747 0.29303279
+[19,] 0.83195027 0.158786065 0.31096864
+[20,] 0.99487195 0.167078366 0.32932628
+[21,] 1.00385049 0.175426033 0.34753107
+[22,] 1.34699176 0.184027027 0.36689077
+[23,] 1.52725592 0.194206812 0.38927550
+[24,] 1.56248132 0.204500082 0.41132235
+[25,] 1.85467950 0.217642436 0.43739727
+[26,] 1.93791303 0.232103772 0.46459689
+[27,] 2.08096688 0.248451634 0.49291950
+[28,] 2.22551104 0.267544962 0.52293238
+[29,] 2.23866481 0.286966924 0.55219966
+[30,] 2.99739926 0.314891008 0.57604023
+[31,] 5.24864987 0.389042862 0.57779641
+attr(,"class")
+[1] "predict.crr"
+> plot(wp)
+> plot(wp,lty=1,col=c(2,4))
+> d2 <- cbind(d2,cv3=cv[,3],cv1=cv[,1])
+> attach(d2)
+> print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cov2,tf=tf,cengroup=cv3))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.02944 -0.06263 -0.42160
+ cov21*uft cov22*tf2
+ -0.84780 0.29850
+standard errors:
+[1] 0.8029 0.3583 0.3595 1.3090 0.4004
+two-sided p-values:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.97 0.86 0.24
+ cov21*uft cov22*tf2
+ 0.52 0.46
+> print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cov2,tf=tf))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.03695 -0.05995 -0.46280
+ cov21*uft cov22*tf2
+ -0.89860 0.32020
+standard errors:
+[1] 0.7862 0.3579 0.3538 1.2180 0.3547
+two-sided p-values:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.96 0.87 0.19
+ cov21*uft cov22*tf2
+ 0.46 0.37
+> print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cbind(cv1,cv1),tf=tf,cengroup=cv3))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.02944 -0.06263 -0.42160
+ cv1*uft cv1*tf2
+ -0.84780 0.29850
+standard errors:
+[1] 0.8029 0.3583 0.3595 1.3090 0.4004
+two-sided p-values:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.97 0.86 0.24
+ cv1*uft cv1*tf2
+ 0.52 0.46
+> print(ww <- crr(ssd,ccd,cbind(cv[,1:2],cv3),cbind(cv1,cv1),tf=tf,cengroup=cv3,subset=X == 1))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.9918 0.1967 -0.2967
+ cv1*uft cv1*tf2
+ -1.8020 0.4035
+standard errors:
+[1] 1.0290 0.4336 0.4546 1.3700 0.3009
+two-sided p-values:
+cbind(cv[, 1:2], cv3)1 cbind(cv[, 1:2], cv3)2 cv3
+ 0.33 0.65 0.51
+ cv1*uft cv1*tf2
+ 0.19 0.18
+> print(summary(ww))
+Competing Risks Regression
+
+Call:
+crr(ftime = ssd, fstatus = ccd, cov1 = cbind(cv[, 1:2], cv3),
+ cov2 = cbind(cv1, cv1), tf = tf, cengroup = cv3, subset = X ==
+ 1)
+
+ coef exp(coef) se(coef) z p-value
+cbind(cv[, 1:2], cv3)1 0.992 2.696 1.029 0.964 0.33
+cbind(cv[, 1:2], cv3)2 0.197 1.217 0.434 0.454 0.65
+cv3 -0.297 0.743 0.455 -0.653 0.51
+cv1*uft -1.802 0.165 1.370 -1.315 0.19
+cv1*tf2 0.403 1.497 0.301 1.341 0.18
+
+ exp(coef) exp(-coef) 2.5% 97.5%
+cbind(cv[, 1:2], cv3)1 2.696 0.371 0.3590 20.25
+cbind(cv[, 1:2], cv3)2 1.217 0.821 0.5204 2.85
+cv3 0.743 1.345 0.3049 1.81
+cv1*uft 0.165 6.064 0.0112 2.42
+cv1*tf2 1.497 0.668 0.8300 2.70
+
+Num. cases = 77 (3 cases omitted due to missing values)
+Pseudo Log-likelihood = -79.2
+Pseudo likelihood ratio test = 2.53 on 5 df,
+> detach(d2)
+> print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3],subset=d2$X==1))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.9918 0.1967 -0.2967 -1.8020 0.4035
+standard errors:
+[1] 1.0290 0.4336 0.4546 1.3700 0.3009
+two-sided p-values:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.33 0.65 0.51 0.19 0.18
+> print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3],failcode=2))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ -0.05102 -0.43670 0.38270 0.93320 -0.41840
+standard errors:
+[1] 0.7200 0.3379 0.3738 1.2430 0.4037
+two-sided p-values:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.94 0.20 0.31 0.45 0.30
+> print(ww <- crr(ss,cc,cv,cov2,tf=tf,cengroup=cv[,3],cencode=2))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.06643 -0.16260 -0.19010 -0.86200 0.32170
+standard errors:
+[1] 0.7965 0.3492 0.3795 1.4070 0.4516
+two-sided p-values:
+ cv1 cv2 cv3 cov21*uft cov22*tf2
+ 0.93 0.64 0.62 0.54 0.48
+> print(ww <- crr(ss,cc,cv[,1]))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+cv[, 1]1
+ -0.1753
+standard errors:
+[1] 0.354
+two-sided p-values:
+cv[, 1]1
+ 0.62
+> print(ww <- crr(ss,cc,cov2=cv[,1],tf=function(x) x))
+3 cases omitted due to missing values
+convergence: TRUE
+coefficients:
+cv[, 1]1*tf1
+ 0.007688
+standard errors:
+[1] 0.2191
+two-sided p-values:
+cv[, 1]1*tf1
+ 0.97
+> print(summary(ww))
+Competing Risks Regression
+
+Call:
+crr(ftime = ss, fstatus = cc, cov2 = cv[, 1], tf = function(x) x)
+
+ coef exp(coef) se(coef) z p-value
+cv[, 1]1*tf1 0.00769 1.01 0.219 0.0351 0.97
+
+ exp(coef) exp(-coef) 2.5% 97.5%
+cv[, 1]1*tf1 1.01 0.992 0.656 1.55
+
+Num. cases = 97 (3 cases omitted due to missing values)
+Pseudo Log-likelihood = -125
+Pseudo likelihood ratio test = 0 on 1 df,
+>
+> proc.time()
+ user system elapsed
+ 0.424 0.026 0.442
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-cmprsk.git
More information about the debian-med-commit
mailing list