R Under development (unstable) (2024-12-12 r87438 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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.

> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
> 
> #
> # Test some more features of surv.diff
> #
> #  First, what happens when one group is a dummy
> #
> 
> 
> #
> # The AML data, with a third group of early censorings "tacked on"
> #
> aml3 <- list(time=   c( 9, 13, 13, 18, 23, 28, 31, 34, 45, 48, 161,
+ 			   5,  5,  8,  8, 12, 16, 23, 27, 30, 33, 43, 45,
+ 			   1,  2,  2,  3,  3,  3,  4),
+ 	    status= c( 1,1,0,1,1,0,1,1,0,1,0, 1,1,1,1,1,0,1,1,1,1,1,1,
+ 				0,0,0,0,0,0,0),
+ 	    x     = as.factor(c(rep("Maintained", 11),
+ 				rep("Nonmaintained", 12), rep("Dummy",7) )))
> 
> aml3 <- data.frame(aml3)
> 
> # These should give the same result (chisq, df), but the second has an
> #  extra group
> survdiff(Surv(time, status) ~x, aml)
Call:
survdiff(formula = Surv(time, status) ~ x, data = aml)

                 N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained    11        7    10.69      1.27       3.4
x=Nonmaintained 12       11     7.31      1.86       3.4

 Chisq= 3.4  on 1 degrees of freedom, p= 0.07 
> survdiff(Surv(time, status) ~x, aml3)
Call:
survdiff(formula = Surv(time, status) ~ x, data = aml3)

                 N Observed Expected (O-E)^2/E (O-E)^2/V
x=Dummy          7        0     0.00       NaN       NaN
x=Maintained    11        7    10.69      1.27       3.4
x=Nonmaintained 12       11     7.31      1.86       3.4

 Chisq= 3.4  on 1 degrees of freedom, p= 0.07 
> 
> 
> #
> # Now a test of the stratified log-rank
> #   There are no tied times within institution, so the coxph program
> #   can be used to give a complete test
> #
> fit <- survdiff(Surv(time, status) ~ pat.karno + strata(inst), lung)
> 
> cfit <- coxph(Surv(time, status) ~ factor(pat.karno) + strata(inst),
+ 		lung, iter=0)
> 
> tdata <- na.omit(lung[,c('time', 'status', 'pat.karno', 'inst')])
> 
> temp1 <- tapply(tdata$status-1, list(tdata$pat.karno, tdata$inst), sum)
> temp1 <- ifelse(is.na(temp1), 0, temp1)
> temp2 <- tapply(cfit$resid,  list(tdata$pat.karno, tdata$inst), sum)
> temp2 <- ifelse(is.na(temp2), 0, temp2)
> 
> temp2 <- temp1 - temp2
> 
> #Now temp1=observed, temp2=expected
> all.equal(c(temp1), c(fit$obs))
[1] TRUE
> all.equal(c(temp2), c(fit$exp))
[1] TRUE
> 
> all.equal(fit$var[-1,-1], solve(cfit$var))
[1] TRUE
> 
> rm(tdata, temp1, temp2)
> 
> proc.time()
   user  system elapsed 
   0.81    0.10    0.85