[R] command in survival package

[Marc Schwartz (via MN)]

On Mon, 2006-01-23 at 10:52 -0800, Thomas Lumley wrote:
> On Mon, 23 Jan 2006, Linda Lei wrote:
> 
> > Thank you guys.
> > But I tried the commands and I still get:
> >
> >> aml1<-aml[aml$group==1,]
> >> aml1
> > [1] time   status x
> > <0 rows> (or 0-length row.names)
> >> esf.fit <- survfit(Surv(aml1$weeks,status) ~ 1)
> > Error in Surv(aml1$weeks, status) : Time variable is not numeric
> > In addition: Warning message:
> > is.na() applied to non-(list or vector) in: is.na(time)
> >
> > which still looks confusing. Or are they should be applied in s-plus
> > instead of R?
> 
> You *still* haven't said where the aml dataset comes from, but as I said 
> before, it isn't the same as the one built in to R, so you need to load it 
> somehow.  The one in R has columns called "time", "status" and "x"; you 
> appear to want one with columns "weeks","status", and "group", and with 
> "group" being numeric rather than factor.
> 
> 
>  	-thomas


I don't have the book with me (it's at home), but I suspect that the
problem here is that the 'aml' dataset provided at the aforementioned
web site: 

  http://www.crcpress.com/e_products/downloads/download.asp?cat_no=C4088

which is to be used with the code under discussion, was more than likely
not imported into an R session by Linda.

Since one needs to do a:

  > library(survival) 

to utilize the R code examples, that loads that package's 'aml' dataset,
resulting in:

  > aml[aml$group==1,]
[1] time   status x     
<0 rows> (or 0-length row.names)

since there is no 'group' column in the built in 'aml' dataset as Thomas
points out.

There is such a structure to the authors' version of the dataset:

> str(read.table("aml.txt", header = TRUE))
`data.frame':	23 obs. of  3 variables:
 $ weeks : num  9 13 13 18 23 28 31 34 45 48 ...
 $ group : num  1 1 1 1 1 1 1 1 1 1 ...
 $ status: num  1 1 0 1 1 0 1 1 0 1 ...


Thus:

>  aml <- read.table("aml.txt", header = TRUE)

> (aml1 <- aml[aml$group==1,])
   weeks group status
1      9     1      1
2     13     1      1
3     13     1      0
4     18     1      1
5     23     1      1
6     28     1      0
7     31     1      1
8     34     1      1
9     45     1      0
10    48     1      1
11   161     1      0


This bring us to the still remaining error in the authors' survival
code, which is not fully corrected in the authors' online errata:

> esf.fit <- survfit(Surv(aml1$weeks,status) ~ 1)
Error in Surv(aml1$weeks, status) : object "status" not found


Thus, something like one of the following should be used instead:

> esf.fit <- survfit(Surv(aml1$weeks, aml1$status) ~ 1)

> esf.fit <- with(aml1, survfit(Surv(weeks, status) ~ 1))


> summary(esf.fit)
Call: survfit(formula = Surv(aml1$weeks, aml1$status) ~ 1)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    9     11       1    0.909  0.0867       0.7541        1.000
   13     10       1    0.818  0.1163       0.6192        1.000
   18      8       1    0.716  0.1397       0.4884        1.000
   23      7       1    0.614  0.1526       0.3769        0.999
   31      5       1    0.491  0.1642       0.2549        0.946
   34      4       1    0.368  0.1627       0.1549        0.875
   48      2       1    0.184  0.1535       0.0359        0.944


I'll pass on a note to the authors, whose e-mail addresses I have
located via Google.

HTH,

Marc Schwartz