---

LEVIS KIBET

21/06461

“ASSIGN 3”

library(survival)
library(survminer)

## Loading required package: ggplot2

## Loading required package: ggpubr

## 
## Attaching package: 'survminer'

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

##data

my_data<-ovarian;my_data

##    futime fustat     age resid.ds rx ecog.ps
## 1      59      1 72.3315        2  1       1
## 2     115      1 74.4932        2  1       1
## 3     156      1 66.4658        2  1       2
## 4     421      0 53.3644        2  2       1
## 5     431      1 50.3397        2  1       1
## 6     448      0 56.4301        1  1       2
## 7     464      1 56.9370        2  2       2
## 8     475      1 59.8548        2  2       2
## 9     477      0 64.1753        2  1       1
## 10    563      1 55.1781        1  2       2
## 11    638      1 56.7562        1  1       2
## 12    744      0 50.1096        1  2       1
## 13    769      0 59.6301        2  2       2
## 14    770      0 57.0521        2  2       1
## 15    803      0 39.2712        1  1       1
## 16    855      0 43.1233        1  1       2
## 17   1040      0 38.8932        2  1       2
## 18   1106      0 44.6000        1  1       1
## 19   1129      0 53.9068        1  2       1
## 20   1206      0 44.2055        2  2       1
## 21   1227      0 59.5890        1  2       2
## 22    268      1 74.5041        2  1       2
## 23    329      1 43.1370        2  1       1
## 24    353      1 63.2192        1  2       2
## 25    365      1 64.4247        2  2       1
## 26    377      0 58.3096        1  2       1

fit kaplan meier

km_fit <- survfit(Surv(futime, fustat) ~ rx, data = my_data)

plot kaplan meier curve

plot(km_fit, xlab = "Time (days)", ylab = "Survival probability", col = c("blue", "red"), lty = 1:2)
legend("topright", legend = c("Treatment 1", "Treatment 2"), col = c("blue", "red"), lty = 1:2)

fit nelson aalen

na_fit <- survfit(Surv(futime, fustat) ~ 1, data =my_data, type="fh")

plot nelson aalen

plot(na_fit, fun = "cumhaz", xlab = "Time (days)", ylab = "Cumulative Hazard", 
     main = "Nelson-Aalen Cumulative Hazard Curve", col = "blue", lwd = 2)

log rank test

log_rank_test <- survdiff(Surv(futime, fustat) ~ rx, data = my_data)

log rank test results

print(log_rank_test)

## Call:
## survdiff(formula = Surv(futime, fustat) ~ rx, data = my_data)
## 
##       N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=1 13        7     5.23     0.596      1.06
## rx=2 13        5     6.77     0.461      1.06
## 
##  Chisq= 1.1  on 1 degrees of freedom, p= 0.3

interpretation of log rank results

The p-value of 0.3 which is greater than 0.05 thus indicates that there is not enough evidence to reject the null hypothesis at a conventional significance level . This means that we do not have sufficient evidence to conclude that there is a significant difference in survival between the two groups

fit cox proportional hazard model

cox_fit <- coxph(Surv(futime, fustat) ~ rx, data = my_data)

cox results

summary(cox_fit)

## Call:
## coxph(formula = Surv(futime, fustat) ~ rx, data = my_data)
## 
##   n= 26, number of events= 12 
## 
##       coef exp(coef) se(coef)      z Pr(>|z|)
## rx -0.5964    0.5508   0.5870 -1.016     0.31
## 
##    exp(coef) exp(-coef) lower .95 upper .95
## rx    0.5508      1.816    0.1743      1.74
## 
## Concordance= 0.608  (se = 0.07 )
## Likelihood ratio test= 1.05  on 1 df,   p=0.3
## Wald test            = 1.03  on 1 df,   p=0.3
## Score (logrank) test = 1.06  on 1 df,   p=0.3

interpretation of cox results

the results suggest that there is no statistically significant association between the treatment or exposure variable and the hazard of the event. The model has a moderate ability to rank the survival times correctly, but the overall fit of the model is not statistically significant.

cox_zph <- cox.zph(cox_fit)
print(cox_zph)

##        chisq df   p
## rx      2.68  1 0.1
## GLOBAL  2.68  1 0.1

plot cox

plot(cox_zph)