R Notebook

library(survival)
library(flexsurv)

## Warning: package 'flexsurv' was built under R version 4.0.3

G1=c(12.3,5.4,12.5,8.6,6.7,9.8,4.3)
G1=sort(G1)
G1

## [1]  4.3  5.4  6.7  8.6  9.8 12.3 12.5

G2=c(2.8,4.3,10.1,5.4,6.7,5.4,9.2,7.2)

G2=sort(G2)
G2

## [1]  2.8  4.3  5.4  5.4  6.7  7.2  9.2 10.1

S1=c(1,1,0,1,1,0,0)
S2=c(1,1,1,1,1,1,0,0)
Surv1=Surv(time=G1,event=S1)
Surv1

## [1]  4.3   5.4   6.7+  8.6   9.8  12.3+ 12.5+

Surv2=Surv(time=G2,event=S2)
Surv2

## [1]  2.8   4.3   5.4   5.4   6.7   7.2   9.2+ 10.1+

s11=survfit(Surv1~1,type="kaplan-meier")
H01=s11$cumhaz
s11$surv

## [1] 0.8571429 0.7142857 0.7142857 0.5357143 0.3571429 0.3571429 0.3571429

H01

## [1] 0.1428571 0.3095238 0.3095238 0.5595238 0.8928571 0.8928571 0.8928571

h01=1:7
h01[1]=H01[1]
for(i in 2:7){
  h01[i]=H01[i]-H01[i-1]
}
h01

## [1] 0.1428571 0.1666667 0.0000000 0.2500000 0.3333333 0.0000000 0.0000000

plot(Surv1,col = "red",ylim=c(0,1.5),xlab="Time",ylab="S(t)",main="Survival Function of Group 1")

plot(G1,h01,type="l",main="Hazard Function of Group 1")

s22=survfit(Surv2~1,type="kaplan-meier")
s22

## Call: survfit(formula = Surv2 ~ 1, type = "kaplan-meier")
## 
##       n  events  median 0.95LCL 0.95UCL 
##    8.00    6.00    6.05    5.40      NA

H02=s22$cumhaz
s22$surv

## [1] 0.875 0.750 0.500 0.375 0.250 0.250 0.250

s22$time

## [1]  2.8  4.3  5.4  6.7  7.2  9.2 10.1

H02

## [1] 0.1250000 0.2678571 0.6011905 0.8511905 1.1845238 1.1845238 1.1845238

h02=1:8
h02[1]=H02[1]
for(i in 2:8){
  h02[i]=H02[i]-H02[i-1]
}
h02

## [1] 0.1250000 0.1428571 0.3333333 0.2500000 0.3333333 0.0000000 0.0000000
## [8]        NA

plot(Surv2,col = "red",ylim=c(0,1.5),xlab="Time",ylab="S(t)",main="Survival Function of Group 1")

plot(G2,h02,type="l",main="Hazard Function of Group 1")

data2=data.frame(cbind(c(G1,G2),c(S1,S2)))
group=c(rep("G1",7),rep("G2",8))
data2["group"]=group
data2

s3=Surv(time=data2$X1,event=data2$X2,type="right")
survdiff(s3~group,data=data2)

## Call:
## survdiff(formula = s3 ~ group, data = data2)
## 
##          N Observed Expected (O-E)^2/E (O-E)^2/V
## group=G1 7        4     5.51     0.414      1.02
## group=G2 8        6     4.49     0.508      1.02
## 
##  Chisq= 1  on 1 degrees of freedom, p= 0.3

data3=data.frame(G1,S1)
Sd1=Surv(G1,S1)
library(frailtypack)

## Warning: package 'frailtypack' was built under R version 4.0.4

## Loading required package: boot

## 
## Attaching package: 'boot'

## The following object is masked from 'package:survival':
## 
##     aml

## Loading required package: MASS

## Loading required package: survC1

## Warning: package 'survC1' was built under R version 4.0.4

## Loading required package: doBy

## Warning: package 'doBy' was built under R version 4.0.4

## 
## Attaching package: 'frailtypack'

## The following object is masked from 'package:survival':
## 
##     cluster

s23=frailtyPenal(Sd1~1,n.knots = 9,kappa = 1,data=data3)

## 
## Be patient. The program is computing ... 
## The program took 0.03 seconds

plot(s23,main="Frailty Hazard Curve of Group 1")

data3=data.frame(G1=2,S2)
Sd1=Surv(G2,S2)
library(frailtypack)
s23=frailtyPenal(Sd1~1,n.knots = 9,kappa = 1,data=data3)

## 
## Be patient. The program is computing ... 
## The program took 0.01 seconds

plot(s23, main="Frailty Hazard Curve of Group 2")