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("g03.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class542")
> setwd("~/Taught1/960-542/Data")
> #***************************************************/
> # Survival function Intervals after Transformation */
> #***************************************************/
> # Employ the log log transformation to keep intervals strictly
> # between 0 and 1.
> library(survival)#For survfit and Surv
> #*******************************************************/
> # 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)
> cat('\n Larynx survival with conf. limits, log log scale  \n')

 Larynx survival with conf. limits, log log scale  
> plot(
+    survfit(Surv(time,delta)~1,conf.type="log-log",data=larynx),
+    main="Larynx Survival", xlab="Time (months)",ylab="Survival",
+    sub="Intervals with log log transform")
> cat('\n Weaning data for mothers who smoked and drank  \n')

 Weaning data for mothers who smoked and drank  
> #***************************************************/
> # Weaning data from Klein and Moeschberger.       */
> # See class web page for link.  Variables are      */
> #duration of breast feeding, weeks                 */
> #Completed breast feeding? (1=yes, 0=no)           */
> #Race of mother (1=white, 2=black, 3=other)        */
> #Mother in poverty (1=yes, 0=no)                   */
> #Mother smoked at birth of child (1=yes, 0=no)     */
> #Mother used alcohol at birth of child (1=yes 0=no)*/
> #Age of mother at birth of child                   */
> #Year of birth                                     */
> #Education level of mother (years of school)       */
> #Prenatal care after 3rd month (1=yes, 0=no)       */
> #***************************************************/
> bfeed<-read.table("bfeed.txt",header=TRUE)
> weanoutl<-survfit(Surv(duration,delta)~1,data=bfeed,
+    subset=(smoke==1)&(alcohol==1), 
+    conf.type="log-log")
> plot(weanoutl,main="Weaning times",
+    xlab="Weaning time (days)",ylab="Survival", sub="log-log ci")
> #Next library must be at least version 1.7; you may need to download 
> #from course web page.
> library(PHInfiniteEstimates)#For simultaneouscoverage
> simultaneouscoverage(1000,20)
[1] 0.747
> #**********************************/
> # Estimators of cumulative hazard.*/
> #**********************************/
> #**************************************/
> # Calculate the Nelson Aalen Estimator*/
> #**************************************/
> library(survival)#For basehaz, quantile.survfit, summary.survfit
> out<-basehaz(coxph(Surv(duration,delta)~1,
+    subset=(smoke==1)&(alcohol==1), data=bfeed))
> #*************************************************/
> # Calculate the log of the Kaplan-Meier Estimator*/
> #*************************************************/
> big<-survfit(Surv(duration,delta)~1,
+    subset=(smoke==1)&(alcohol==1), data=bfeed)
> big$nls<--log(big$surv)
> #*************************************************/
> # Plot both                                      */
> #*************************************************/
> plot(out$time,out$hazard,pch=1,xlab="Time",
+   ylab="Cumulative Hazard")
> points( big$time, big$nls ,pch=2)
> legend(mean(big$time),median(big$nls),pch=1:2,
+   legend=c("Nelson-Allen","log of Kaplan-Meier"))
> #***********************************************/
> # Quantile estimation                          */
> #***********************************************/
> cat('\n Quantile CIs for weaning, high risk  \n')

 Quantile CIs for weaning, high risk  
> library(survival)#For survfit
> #To get other quantiles from quantile.survfit  
> #use argument probs.
> #***************************************************/
> # Weaning data from Klein and Moeschberger.       */
> # See class web page for link.  Variables are      */
> #duration of breast feeding, weeks                 */
> #Completed breast feeding? (1=yes, 0=no)           */
> #Race of mother (1=white, 2=black, 3=other)        */
> #Mother in poverty (1=yes, 0=no)                   */
> #Mother smoked at birth of child (1=yes, 0=no)     */
> #Mother used alcohol at birth of child (1=yes 0=no)*/
> #Age of mother at birth of child                   */
> #Year of birth                                     */
> #Education level of mother (years of school)       */
> #Prenatal care after 3rd month (1=yes, 0=no)       */
> #***************************************************/
> bfeed<-read.table("bfeed.txt",header=TRUE)
> plot(fit1<-survfit(Surv(duration,delta)~1,
+    subset=(alcohol>0)&(smoke>0), data=bfeed), 
+    main="Weaning Data", sub="High Risk Subgroup")
> abline(h=0.5,lty=2)
> # Putting the next line in parentheses triggers printing
> # of computed object.
> (qout<-quantile(fit1,data=bfeed[(alcohol>0)&(smoke>0),]))
$quantile
  25   50   75 
 2.5  5.5 16.0 

$lower
25 50 75 
 1  4 10 

$upper
25 50 75 
 4 12 32 

> abline(v=qout$lower["50"],lty=3)
> abline(v=qout$upper["50"],lty=3)
> #*******************************************************/
> # 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)
> stage1<-larynx[larynx$stage==1,]
> cat('\n Quantile CIs for larynx cancer, low risk  \n')

 Quantile CIs for larynx cancer, low risk  
> #calls quantile.survfit.  To get other quantiles use
> #argument probs.                                    
> fit2<-survfit(Surv(time,delta)~1,data=stage1)
> plot(fit2,main="Larynx Cancer Data",sub="Low Risk")
> abline(h=0.5,lty=2)
> quantile(fit2)
$quantile
 25  50  75 
4.0 6.5  NA 

$lower
 25  50  75 
3.3 5.3 7.4 

$upper
 25  50  75 
6.5  NA  NA 

> #***********************************************/
> # Estimate expectation of lifetime.            */
> #***********************************************/
> #Estimate expectation by stage.  Specify grouping using
> #formula with a numerical value.  This has the effect of
> #stratification, even if strata() is not used.  Summary 
> #gets evaluated using summary.survfit.
> #Package survival doesn't make it easy to get mean estimates
> #out.  Because in many cases the expectation from the KM
> #survival estimate is infinite, R insists on restricted mean.
> #rmean controls this.  Set it to a numerical value to set
> #the restriction.  Set it to "none" (the default) to suppress
> #calculation.  Setting it to rmean="common" truncates survival 
> #curves not finishing at 0 to the maximum time represented 
> #in the data set. Setting it to "individual" sets a separate 
> #restricted mean for each group.  I don't see why.  If rmean
> #is calculated, it is column 5 of the resulting table.
> summary(survfit(Surv(time,delta)~stage,data=larynx),
+     rmean="common")$table[,5]
 stage=1  stage=2  stage=3  stage=4 
7.083285 6.866786 5.127222 2.564103 
>