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("g06.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class542")
> setwd("~/Taught1/960-542/Data")
> library(survival)#For coxph
> #*******************************************************/
> # 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)
> #**************************************************************/
> # Fit proportional hazards model.  Stage gets converted to a  */ 
> # factor before removing stages 1 and 4.  Stage 1 gets treated*/
> # as baseline, even though it's empty, and stage 4 gets       */
> # estimated as NA since it's empty.                           */
> #**************************************************************/
> summary(coxph(Surv(time,delta)~factor(stage),
+   data=larynx,subset=(stage==2)|(stage==3),method="breslow"))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    subset = (stage == 2) | (stage == 3), method = "breslow")

  n= 44, number of events= 24 

                  coef exp(coef) se(coef)      z Pr(>|z|)
factor(stage)2 -0.5386    0.5835   0.4514 -1.193    0.233
factor(stage)3      NA        NA   0.0000     NA       NA
factor(stage)4      NA        NA   0.0000     NA       NA

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2    0.5835      1.714    0.2409     1.413
factor(stage)3        NA         NA        NA        NA
factor(stage)4        NA         NA        NA        NA

Concordance= 0.585  (se = 0.051 )
Likelihood ratio test= 1.52  on 1 df,   p=0.2
Wald test            = 1.42  on 1 df,   p=0.2
Score (logrank) test = 1.46  on 1 df,   p=0.2

> #***********************************************************/
> # Policies for ties.  Compare earlier lifetest results. Log*/
> # rank statistic corresponds to score statistic for reg-   */
> # ression model.                                           */
> # Default tie handling in SAS is Breslow and in R is Efron.*/
> #***********************************************************/
> #R survdiff has no tie option.  
> cat(tl<-"Larynx Cancer Survival Stages 2&3\n",
+   "Testing Stage")
Larynx Cancer Survival Stages 2&3
 Testing Stage> survdiff(Surv(time,delta)~factor(stage),data=larynx,
+    subset=(stage==2)|(stage==3))
Call:
survdiff(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    subset = (stage == 2) | (stage == 3))

                 N Observed Expected (O-E)^2/E (O-E)^2/V
factor(stage)=2 17        7     9.89     0.846      1.47
factor(stage)=3 27       17    14.11     0.593      1.47

 Chisq= 1.5  on 1 degrees of freedom, p= 0.2 
> #*******************************************************/
> #Breslow option                                        */
> #*******************************************************/
> cat(tl,"Tie handling Breslow\n")
Larynx Cancer Survival Stages 2&3
 Tie handling Breslow
> summary(coxph(Surv(time,delta)~factor(stage),data=larynx,
+    method="breslow", subset=(stage==2)|(stage==3)))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    subset = (stage == 2) | (stage == 3), method = "breslow")

  n= 44, number of events= 24 

                  coef exp(coef) se(coef)      z Pr(>|z|)
factor(stage)2 -0.5386    0.5835   0.4514 -1.193    0.233
factor(stage)3      NA        NA   0.0000     NA       NA
factor(stage)4      NA        NA   0.0000     NA       NA

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2    0.5835      1.714    0.2409     1.413
factor(stage)3        NA         NA        NA        NA
factor(stage)4        NA         NA        NA        NA

Concordance= 0.585  (se = 0.051 )
Likelihood ratio test= 1.52  on 1 df,   p=0.2
Wald test            = 1.42  on 1 df,   p=0.2
Score (logrank) test = 1.46  on 1 df,   p=0.2

> #*******************************************************/
> #Efron option                                          */
> #*******************************************************/
> cat(tl,"Ties handled using Efron method\n")
Larynx Cancer Survival Stages 2&3
 Ties handled using Efron method
> summary(coxph(Surv(time,delta)~factor(stage),
+   data=larynx,method="efron",subset=(stage==2)|(stage==3)))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    subset = (stage == 2) | (stage == 3), method = "efron")

  n= 44, number of events= 24 

                  coef exp(coef) se(coef)      z Pr(>|z|)
factor(stage)2 -0.5425    0.5813   0.4514 -1.202    0.229
factor(stage)3      NA        NA   0.0000     NA       NA
factor(stage)4      NA        NA   0.0000     NA       NA

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2    0.5813       1.72      0.24     1.408
factor(stage)3        NA         NA        NA        NA
factor(stage)4        NA         NA        NA        NA

Concordance= 0.585  (se = 0.051 )
Likelihood ratio test= 1.54  on 1 df,   p=0.2
Wald test            = 1.44  on 1 df,   p=0.2
Score (logrank) test = 1.48  on 1 df,   p=0.2

> # Efron is the default value.
> summary(coxph(Surv(time,delta)~factor(stage),
+   data=larynx,subset=(stage==2)|(stage==3)))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    subset = (stage == 2) | (stage == 3))

  n= 44, number of events= 24 

                  coef exp(coef) se(coef)      z Pr(>|z|)
factor(stage)2 -0.5425    0.5813   0.4514 -1.202    0.229
factor(stage)3      NA        NA   0.0000     NA       NA
factor(stage)4      NA        NA   0.0000     NA       NA

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2    0.5813       1.72      0.24     1.408
factor(stage)3        NA         NA        NA        NA
factor(stage)4        NA         NA        NA        NA

Concordance= 0.585  (se = 0.051 )
Likelihood ratio test= 1.54  on 1 df,   p=0.2
Wald test            = 1.44  on 1 df,   p=0.2
Score (logrank) test = 1.48  on 1 df,   p=0.2

> #*******************************************************/
> # Discrete/ Exact option.                              */
> #*******************************************************/
> cat(tl,"Tie handling Exact\n")
Larynx Cancer Survival Stages 2&3
 Tie handling Exact
> summary(coxph(Surv(time,delta)~factor(stage),
+   data=larynx,method="exact",subset=(stage==2)|(stage==3)))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    subset = (stage == 2) | (stage == 3), method = "exact")

  n= 44, number of events= 24 

                  coef exp(coef) se(coef)      z Pr(>|z|)
factor(stage)2 -0.5424    0.5814   0.4529 -1.198    0.231
factor(stage)3      NA        NA   0.0000     NA       NA
factor(stage)4      NA        NA   0.0000     NA       NA

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2    0.5814       1.72    0.2393     1.412
factor(stage)3        NA         NA        NA        NA
factor(stage)4        NA         NA        NA        NA

Concordance= 0.585  (se = 0.051 )
Likelihood ratio test= 1.53  on 1 df,   p=0.2
Wald test            = 1.43  on 1 df,   p=0.2
Score (logrank) test = 1.47  on 1 df,   p=0.2

> #************************************************/
> #K sample testing via dummy variables.          */
> #************************************************/
> cat(tl)
Larynx Cancer Survival Stages 2&3
> summary(coxph(Surv(time,delta)~factor(stage),
+   data=larynx,method="breslow"))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx, 
    method = "breslow")

  n= 90, number of events= 50 

                  coef exp(coef) se(coef)     z Pr(>|z|)    
factor(stage)2 0.06576   1.06797  0.45844 0.143   0.8859    
factor(stage)3 0.61206   1.84423  0.35520 1.723   0.0849 .  
factor(stage)4 1.72284   5.60040  0.41966 4.105 4.04e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2     1.068     0.9364    0.4348     2.623
factor(stage)3     1.844     0.5422    0.9193     3.700
factor(stage)4     5.600     0.1786    2.4604    12.748

Concordance= 0.668  (se = 0.037 )
Likelihood ratio test= 16.26  on 3 df,   p=0.001
Wald test            = 18.95  on 3 df,   p=3e-04
Score (logrank) test = 22.46  on 3 df,   p=5e-05

> # By default in R reference category is the first.
> cat("Change baseline in coxph\n")
Change baseline in coxph
> summary(coxph(Surv(time,delta)~factor(stage,
+   levels=c(4,3,2,1)), data=larynx,method="breslow"))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage, levels = c(4, 
    3, 2, 1)), data = larynx, method = "breslow")

  n= 90, number of events= 50 

                                          coef exp(coef)
factor(stage, levels = c(4, 3, 2, 1))3 -1.1108    0.3293
factor(stage, levels = c(4, 3, 2, 1))2 -1.6571    0.1907
factor(stage, levels = c(4, 3, 2, 1))1 -1.7228    0.1786
                                       se(coef)      z Pr(>|z|)
factor(stage, levels = c(4, 3, 2, 1))3   0.4063 -2.734  0.00625
factor(stage, levels = c(4, 3, 2, 1))2   0.4981 -3.326  0.00088
factor(stage, levels = c(4, 3, 2, 1))1   0.4197 -4.105 4.04e-05
                                          
factor(stage, levels = c(4, 3, 2, 1))3 ** 
factor(stage, levels = c(4, 3, 2, 1))2 ***
factor(stage, levels = c(4, 3, 2, 1))1 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                       exp(coef) exp(-coef)
factor(stage, levels = c(4, 3, 2, 1))3    0.3293      3.037
factor(stage, levels = c(4, 3, 2, 1))2    0.1907      5.244
factor(stage, levels = c(4, 3, 2, 1))1    0.1786      5.600
                                       lower .95 upper .95
factor(stage, levels = c(4, 3, 2, 1))3   0.14852    0.7301
factor(stage, levels = c(4, 3, 2, 1))2   0.07183    0.5063
factor(stage, levels = c(4, 3, 2, 1))1   0.07845    0.4064

Concordance= 0.668  (se = 0.037 )
Likelihood ratio test= 16.26  on 3 df,   p=0.001
Wald test            = 18.95  on 3 df,   p=3e-04
Score (logrank) test = 22.46  on 3 df,   p=5e-05

> #************************************************/
> #Models and Submodels                           */
> #************************************************/
> summary(f1<-coxph(Surv(time,delta)~factor(stage),
+   data=larynx))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx)

  n= 90, number of events= 50 

                  coef exp(coef) se(coef)     z Pr(>|z|)    
factor(stage)2 0.06481   1.06696  0.45843 0.141   0.8876    
factor(stage)3 0.61481   1.84930  0.35519 1.731   0.0835 .  
factor(stage)4 1.73490   5.66838  0.41939 4.137 3.52e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

               exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2     1.067     0.9372    0.4344      2.62
factor(stage)3     1.849     0.5407    0.9219      3.71
factor(stage)4     5.668     0.1764    2.4916     12.90

Concordance= 0.668  (se = 0.037 )
Likelihood ratio test= 16.49  on 3 df,   p=9e-04
Wald test            = 19.24  on 3 df,   p=2e-04
Score (logrank) test = 22.88  on 3 df,   p=4e-05

> cat(tl,"Using stage as quant. instead of categorical\n")
Larynx Cancer Survival Stages 2&3
 Using stage as quant. instead of categorical
> summary(f2<-coxph(Surv(time,delta)~stage,data=larynx))
Call:
coxph(formula = Surv(time, delta) ~ stage, data = larynx)

  n= 90, number of events= 50 

        coef exp(coef) se(coef)     z Pr(>|z|)    
stage 0.5088    1.6633   0.1412 3.604 0.000313 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

      exp(coef) exp(-coef) lower .95 upper .95
stage     1.663     0.6012     1.261     2.193

Concordance= 0.668  (se = 0.037 )
Likelihood ratio test= 13.26  on 1 df,   p=3e-04
Wald test            = 12.99  on 1 df,   p=3e-04
Score (logrank) test = 13.82  on 1 df,   p=2e-04

> cat("Test stage as cts vs. stage as factor\n")
Test stage as cts vs. stage as factor
> #Test adequacy of continuous model via likelihood
> #ratio.
> 1-pchisq(2*(f1$loglik-f2$loglik)[2],2)
[1] 0.199715
> #Compare via the built-in test.
> anova(f2,f1)
Analysis of Deviance Table
 Cox model: response is  Surv(time, delta)
 Model 1: ~ stage
 Model 2: ~ factor(stage)
   loglik  Chisq Df Pr(>|Chi|)
1 -190.23                     
2 -188.62 3.2217  2     0.1997
> #****************************************/
> #Add interactions between stage and age.*/
> #****************************************/
> cat("Age added as dichotomous covariate.\n")
Age added as dichotomous covariate.
> cat(tl,"Adding interactions between age and stage.\n")
Larynx Cancer Survival Stages 2&3
 Adding interactions between age and stage.
> #People in later stages tend to be older.  As with other
> #regression models, this undermines significance of each.
> summary(f4<-coxph(Surv(time,delta)~factor(stage)+factor(age>65),
+    data=larynx))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage) + factor(age > 
    65), data = larynx)

  n= 90, number of events= 50 

                        coef exp(coef) se(coef)     z Pr(>|z|)
factor(stage)2       0.08082   1.08418  0.46351 0.174   0.8616
factor(stage)3       0.62024   1.85937  0.35597 1.742   0.0814
factor(stage)4       1.72984   5.63977  0.42026 4.116 3.85e-05
factor(age > 65)TRUE 0.06824   1.07062  0.28938 0.236   0.8136
                        
factor(stage)2          
factor(stage)3       .  
factor(stage)4       ***
factor(age > 65)TRUE    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                     exp(coef) exp(-coef) lower .95 upper .95
factor(stage)2           1.084     0.9224    0.4371     2.689
factor(stage)3           1.859     0.5378    0.9255     3.736
factor(stage)4           5.640     0.1773    2.4748    12.852
factor(age > 65)TRUE     1.071     0.9340    0.6072     1.888

Concordance= 0.661  (se = 0.039 )
Likelihood ratio test= 16.54  on 4 df,   p=0.002
Wald test            = 19.29  on 4 df,   p=7e-04
Score (logrank) test = 22.93  on 4 df,   p=1e-04

> summary(f3<-coxph(Surv(time,delta)~factor(stage)*factor(age>65),
+    data=larynx))
Call:
coxph(formula = Surv(time, delta) ~ factor(stage) * factor(age > 
    65), data = larynx)

  n= 90, number of events= 50 

                                       coef exp(coef) se(coef)
factor(stage)2                      -0.4771    0.6206   0.6935
factor(stage)3                       0.5686    1.7657   0.5047
factor(stage)4                       1.6380    5.1448   0.6358
factor(age > 65)TRUE                -0.1725    0.8416   0.5214
factor(stage)2:factor(age > 65)TRUE  1.1749    3.2378   0.9416
factor(stage)3:factor(age > 65)TRUE  0.0736    1.0764   0.7127
factor(stage)4:factor(age > 65)TRUE  0.2226    1.2494   0.8189
                                         z Pr(>|z|)   
factor(stage)2                      -0.688  0.49151   
factor(stage)3                       1.127  0.25990   
factor(stage)4                       2.576  0.00999 **
factor(age > 65)TRUE                -0.331  0.74080   
factor(stage)2:factor(age > 65)TRUE  1.248  0.21214   
factor(stage)3:factor(age > 65)TRUE  0.103  0.91776   
factor(stage)4:factor(age > 65)TRUE  0.272  0.78572   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                    exp(coef) exp(-coef)
factor(stage)2                         0.6206     1.6113
factor(stage)3                         1.7657     0.5663
factor(stage)4                         5.1448     0.1944
factor(age > 65)TRUE                   0.8416     1.1882
factor(stage)2:factor(age > 65)TRUE    3.2378     0.3089
factor(stage)3:factor(age > 65)TRUE    1.0764     0.9290
factor(stage)4:factor(age > 65)TRUE    1.2494     0.8004
                                    lower .95 upper .95
factor(stage)2                         0.1594     2.416
factor(stage)3                         0.6567     4.748
factor(stage)4                         1.4797    17.887
factor(age > 65)TRUE                   0.3029     2.338
factor(stage)2:factor(age > 65)TRUE    0.5113    20.501
factor(stage)3:factor(age > 65)TRUE    0.2662     4.352
factor(stage)4:factor(age > 65)TRUE    0.2510     6.220

Concordance= 0.68  (se = 0.038 )
Likelihood ratio test= 18.3  on 7 df,   p=0.01
Wald test            = 20.3  on 7 df,   p=0.005
Score (logrank) test = 24.26  on 7 df,   p=0.001

> cat("Test for interactions\n")
Test for interactions
> anova(f4,f3)
Analysis of Deviance Table
 Cox model: response is  Surv(time, delta)
 Model 1: ~ factor(stage) + factor(age > 65)
 Model 2: ~ factor(stage) * factor(age > 65)
   loglik  Chisq Df Pr(>|Chi|)
1 -188.59                     
2 -187.71 1.7634  3     0.6229
> # How does the p-value behave?
> library(PHInfiniteEstimates)#For testcox
> testcox()
ProgressProgressProgressProgressProgressProgressProgressProgressProgressProgress[1] 0.935
>