R version 4.2.2 Patched (2022-11-10 r83330) -- "Innocent and Trusting"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (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.

  Natural language support but running in an English locale

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.

> pdf("g05.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class542")
> setwd("~/Taught1/960-542/Data")
> #************************************************/
> # Reyni    Testing                              */
> #************************************************/
> #*******************************************************/
> # Larynx cancer data from Klelin and Moeschberger. See */
> # class web page for link.  Cancers are graded on an   */
> # integer scale, with increasing numbers indicating    */
> # more sever.  Variables are stage, time to death, year*/
> # of diagnosis, and status.                            */
> #*******************************************************/
> larynx<-read.table("larynx.txt",header=TRUE)
> larynx$fs<- ((larynx$delta==0)&(larynx$time<4.4))|
+   ((larynx$delta==0)&(larynx$time>7.2))|
+   ((larynx$delta==1)&(larynx$time<1.0))|
+   ((larynx$delta==1)&(larynx$time>4.0))
> #*******************************************************/
> #Two-group examples from the last lecture.             */
> #*******************************************************/
> library(PHInfiniteEstimates)#For gehan.wilcoxon.test
> cat('\n Larynx Cancer Surv. Stage 1,4 Tumors by Stage  \n')

 Larynx Cancer Surv. Stage 1,4 Tumors by Stage  
> plot(survfit(Surv(time,delta)~stage,data=larynx,
+    subset=(stage==1)|(stage==4)), conf.int=T,lty=1:2,
+   main="Larynx Cancer Surv. Stage 1,4 Tumors by Stage")
> gehan.wilcoxon.test(Surv(time,delta)~stage,
+    data=larynx[(larynx$stage==1)|(larynx$stage==4),],
+    plot=TRUE,gehan=FALSE)

	log rank

data:  
= 23.4, p-value = 1.315e-06
alternative hypothesis: two-sided

> cat('\n Larynx Cancer Surv. Stage 2,3 Tumors by Stage  \n')

 Larynx Cancer Surv. Stage 2,3 Tumors by Stage  
> plot(survfit(Surv(time,delta)~stage,data=larynx,
+    subset=(stage==2)|(stage==3)), conf.int=T,lty=1:2,
+    main="Larynx Cancer Surv. Stage 2,3 Tumors by Stage")
> gehan.wilcoxon.test(Surv(time,delta)~stage,
+    data=larynx[(larynx$stage==2)|(larynx$stage==3),],
+    plot=TRUE,gehan=FALSE)

	log rank

data:  
= 1.4683, p-value = 0.2256
alternative hypothesis: two-sided

> cat('\n Larynx Cancer Surv. Crossing Survival Curves  \n')

 Larynx Cancer Surv. Crossing Survival Curves  
> plot(survfit(Surv(time,delta)~fs,data=larynx,
+    subset=(stage==1)|(stage==4)), conf.int=T,lty=1:2,
+    main="Larynx Cancer Surv. Stage 1,4 Tumors by Stage",
+    sub="Fake variable")
> gehan.wilcoxon.test(Surv(time,delta)~fs,
+    data=larynx[(larynx$stage==1)|(larynx$stage==4),],
+    plot=TRUE,gehan=FALSE)

	log rank

data:  
= 0.20298, p-value = 0.6523
alternative hypothesis: two-sided

> #######################################################
> # Fake data to show different ways things can happen  #
> #######################################################
> par(mfrow=c(2,2))
> source("fakereyni.R")#Three large simulated data sets.
> plot(survfit(Surv(time,delta)~grp,data=hugefake0),col=1:2)
> plot(survfit(Surv(time,delta)~grp,data=hugefake1),col=1:2)
> plot(survfit(Surv(time,delta)~grp,data=hugefake2),col=1:2)
> par(mfrow=c(2,2))
> # Example of non-significant log rank and Reyni tests.
> gehan.wilcoxon.test(Surv(time,delta)~grp,data=hugefake0,
+    plot=TRUE,gehan=FALSE)

	log rank

data:  
= 0.324, p-value = 0.5692
alternative hypothesis: two-sided

> # Example of significant log rank test.
> gehan.wilcoxon.test(Surv(time,delta)~grp,data=hugefake1,
+    plot=TRUE,gehan=FALSE)

	log rank

data:  
= 5.4732, p-value = 0.01931
alternative hypothesis: two-sided

> # Example of crossing Reyni tests.
> gehan.wilcoxon.test(Surv(time,delta)~grp,data=hugefake2,
+    plot=TRUE,gehan=FALSE)

	log rank

data:  
= 1.0424, p-value = 0.3073
alternative hypothesis: two-sided

> # Compare log rank, KS, and Renyi tests for Larynx data.
> a1<-survdiff(Surv(time,delta)~stage,data=larynx,
+    subset=(stage==2)|(stage==3))
> stages23<-larynx[(larynx$stage==2)|(larynx$stage==3),]
> #The Kolmogorov-Smirnov tests were formerly in surv2sample. The
> #package is no longer active on CRAN, and I can't make it work.  
> #To install an archived package, one could use devtools, via
> # install.packages(devtools);library(devtools)
> # install_version("surv2sample",version="0.1-2")
> # library(surv2sample)#For surv2.ks
> #a2<-surv2.ks(Surv(stages23$time,stages23$delta),
> #  stages23$stage-1)
> # I couldn't make this work.  Here's a replacement:
> library(MultNonParam)#For twosamplesurvpvs
> #P-values are calculated using permutation.  Recall that this
> #requires the assumption of equality of censoring.
> a2<-twosamplesurvpvs(stages23$time,stages23$delta,
+    stages23$stage,plotme=FALSE)[[1]][1]
> #To get Reyni version of log rank test:
> library(survMisc)#For Reyni test, via ten and comp
> temp<-ten(Surv(time,delta)~stage,data=stages23)
> comp(temp);a3<-attributes(temp)$sup[[5]][1]
                    Q        Var       Z pNorm
1             2.89227    5.71165 1.21020     5
n           112.00000 5628.55357 1.49286     1
sqrtN        17.97103  166.13080 1.39427     4
S1            2.44131    3.00347 1.40867     3
S2            2.37871    2.81093 1.41879     2
FH_p=1_q=1    0.31721    0.19142 0.72502     6
              maxAbsZ        Var      Q pSupBr
1             3.78788    5.71165 1.5849      5
n           135.00000 5628.55357 1.7994      1
sqrtN        22.55941  166.13080 1.7503      2
S1            3.01009    3.00347 1.7369      4
S2            2.92756    2.81093 1.7461      3
FH_p=1_q=1    0.52658    0.19142 1.2036      6
> cat('\n Larynx Cancer Surv. Stage 1,4 Tumors by Stage  \n')

 Larynx Cancer Surv. Stage 1,4 Tumors by Stage  
> plot(survfit(Surv(time,delta)~strata(stage),data=larynx,
+    subset=(stage==1)|(stage==4)), conf.int=T,lty=1:2,
+    main="Larynx Cancer Surv. Stage 1,4 Tumors by Stage")
> b1<-survdiff(Surv(time,delta)~stage,data=larynx,
+    subset=(stage==1)|(stage==4))
> stages14<-larynx[(larynx$stage==1)|(larynx$stage==4),]
> b2<-twosamplesurvpvs(stages14$time,stages14$delta,
+    stages14$stage,plotme=FALSE)[[1]][1]
> temp<-ten(Surv(time,delta)~stage,data=stages14)
> comp(temp);b3<-attributes(temp)$sup[[5]][1]
                    Q        Var      Z pNorm
1          7.7807e+00 2.5871e+00 4.8374     4
n          2.9800e+02 3.7968e+03 4.8362     5
sqrtN      4.7942e+01 9.7506e+01 4.8551     1
S1         6.4096e+00 1.7513e+00 4.8435     2
S2         6.2453e+00 1.6638e+00 4.8418     3
FH_p=1_q=1 9.0081e-01 5.2942e-02 3.9150     6
              maxAbsZ        Var      Q pSupBr
1          7.9095e+00 2.5871e+00 4.9175      3
n          3.0100e+02 3.7968e+03 4.8849      6
sqrtN      4.8564e+01 9.7506e+01 4.9181      2
S1         6.4799e+00 1.7513e+00 4.8966      4
S2         6.3127e+00 1.6638e+00 4.8940      5
FH_p=1_q=1 9.3219e-01 5.2942e-02 4.0514      1
> cat('\n Larynx Cancer Surv. by artificial variable  \n') 

 Larynx Cancer Surv. by artificial variable  
> plot(survfit(Surv(time,delta)~strata(fs),stages14),
+    conf.int=T,lty=1:2,
+    main="Larynx Cancer Surv. Stage 1,4 Tumors by Stage")
> c1<-survdiff(Surv(time,delta)~fs,data=stages14)
> c2<-twosamplesurvpvs(stages14$time,stages14$delta,
+    stages14$fs,plotme=FALSE)[[1]][1]
> temp<-ten(Surv(time,delta)~fs,data=stages14)
> comp(temp);c3<-attributes(temp)$sup[[5]][1]
                    Q        Var        Z pNorm
1            -1.10446    6.00952 -0.45054     2
n           -22.00000 6607.18287 -0.27065     4
sqrtN        -5.71712  189.23934 -0.41560     3
S1           -0.44606    3.22292 -0.24847     5
S2           -0.43224    3.03520 -0.24810     6
FH_p=1_q=1   -0.58917    0.20317 -1.30709     1
              maxAbsZ        Var      Q pSupBr
1             3.30152    6.00952 1.3468      6
n           144.00000 6607.18287 1.7716      2
sqrtN        21.80098  189.23934 1.5848      5
S1            3.03951    3.22292 1.6931      4
S2            2.97141    3.03520 1.7056      3
FH_p=1_q=1    1.04254    0.20317 2.3129      1
> out<-cbind(pchisq(c(a1$chisq,b1$chisq,c1$chisq),1,lower.tail=F),
+ c(a2,b2,c2),c(a3,b3,c3))
> dimnames(out)<-list(c("Stage 2 vs 3","Stage 1 vs 4",
+   "fake"),c("Log rank","Kolmogorov-Smirnov","Reni"))
> print(out)
                 Log rank Kolmogorov-Smirnov         Reni
Stage 2 vs 3 2.256073e-01             0.3944 2.259523e-01
Stage 1 vs 4 1.315474e-06             0.0007 1.753388e-06
fake         6.523236e-01             0.2650 3.559995e-01
> #************************************************/
> # Power calculations for the log rank statistic */
> #************************************************/
> library(Hmisc)#For cpower
> #Power of log rank test at 10 years for a study with 100
> #patients per arm, at which half of the control subjects
> #have the event, and 25% fewer treatment patients have
> #the event.  The 0 indicates that subjects all enter the
> #study at the same time (and hence are all followed for
> #the minimum duration), and 10
> # indicates that people are followed for a minimum of 10
> # years.  Exponential survival curves are presumed.
> cpower(10,100,.50,25,0,10)

Accrual duration: 0 y  Minimum follow-up: 10 y

Total sample size: 100 

Alpha= 0.05 

10-year Mortalities
     Control Intervention 
       0.500        0.375 

Hazard Rates
     Control Intervention 
  0.06931472   0.04700036 

Probabilities of an Event During Study
     Control Intervention 
       0.500        0.375 

Expected Number of Events
     Control Intervention 
        25.0         18.8 

Hazard ratio: 0.6780719 
Standard deviation of log hazard ratio: 0.305505 
    Power 
0.2462496 
>