library(survival)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
##
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
##
## myeloma
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
hmohiv<-read.table("https://stats.idre.ucla.edu/stat/r/examples/asa/hmohiv.csv", sep=",", header = TRUE)
hmohiv$s_obj <- Surv(hmohiv$time, hmohiv$censor)
20-29 : cage_1
30-39 : cage_2
40-49 : cage_3
50-54 : cage_4
hmohiv$cage <- car::recode(hmohiv$age, "20:29='cage_1'; 30:39='cage_2'; 40:49='cage_3'; 50:59='cage_4'", as.factor = TRUE)
hmohiv
## ID time age drug censor entdate enddate s_obj cage
## 1 1 5 46 0 1 5/15/1990 10/14/1990 5 cage_3
## 2 2 6 35 1 0 9/19/1989 3/20/1990 6+ cage_2
## 3 3 8 30 1 1 4/21/1991 12/20/1991 8 cage_2
## 4 4 3 30 1 1 1/3/1991 4/4/1991 3 cage_2
## 5 5 22 36 0 1 9/18/1989 7/19/1991 22 cage_2
## 6 6 1 32 1 0 3/18/1991 4/17/1991 1+ cage_2
## 7 7 7 36 1 1 11/11/1989 6/11/1990 7 cage_2
## 8 8 9 31 1 1 11/25/1989 8/25/1990 9 cage_2
## 9 9 3 48 0 1 2/11/1991 5/13/1991 3 cage_3
## 10 10 12 47 0 1 8/11/1989 8/11/1990 12 cage_3
## 11 11 2 28 1 0 4/11/1990 6/10/1990 2+ cage_1
## 12 12 12 34 0 1 5/11/1991 5/10/1992 12 cage_2
## 13 13 1 44 1 1 1/17/1989 2/16/1989 1 cage_3
## 14 14 15 32 1 1 2/16/1991 5/17/1992 15 cage_2
## 15 15 34 36 0 1 4/9/1991 2/6/1994 34 cage_2
## 16 16 1 36 0 1 3/9/1991 4/8/1991 1 cage_2
## 17 17 4 54 0 1 8/3/1990 12/2/1990 4 cage_4
## 18 18 19 35 0 0 6/10/1990 1/8/1992 19+ cage_2
## 19 19 3 44 1 0 6/12/1991 9/11/1991 3+ cage_3
## 20 20 2 38 0 1 1/7/1991 3/8/1991 2 cage_2
## 21 21 2 40 0 0 8/29/1989 10/28/1989 2+ cage_3
## 22 22 6 34 1 1 5/29/1989 11/27/1989 6 cage_2
## 23 23 60 25 0 0 11/16/1990 11/14/1995 60+ cage_1
## 24 24 11 32 0 1 5/9/1990 4/8/1991 11 cage_2
## 25 25 2 42 1 0 9/10/1991 11/9/1991 2+ cage_3
## 26 26 5 47 0 1 12/26/1991 5/26/1992 5 cage_3
## 27 27 4 30 0 0 5/29/1991 9/27/1991 4+ cage_2
## 28 28 1 47 1 1 5/1/1990 5/31/1990 1 cage_3
## 29 29 13 41 0 1 3/24/1991 4/22/1992 13 cage_3
## 30 30 3 40 1 1 7/18/1989 10/17/1989 3 cage_3
## 31 31 2 43 0 1 9/16/1990 11/15/1990 2 cage_3
## 32 32 1 41 0 1 6/22/1989 7/22/1989 1 cage_3
## 33 33 30 30 0 1 4/27/1990 10/25/1992 30 cage_2
## 34 34 7 37 0 1 5/16/1990 12/14/1990 7 cage_2
## 35 35 4 42 1 1 2/19/1989 6/20/1989 4 cage_3
## 36 36 8 31 1 1 2/17/1990 10/18/1990 8 cage_2
## 37 37 5 39 1 1 8/6/1991 1/5/1992 5 cage_2
## 38 38 10 32 0 1 8/10/1989 6/10/1990 10 cage_2
## 39 39 2 51 0 1 12/27/1990 2/25/1991 2 cage_4
## 40 40 9 36 0 1 4/26/1989 1/24/1990 9 cage_2
## 41 41 36 43 0 1 12/4/1990 12/3/1993 36 cage_3
## 42 42 3 39 0 1 4/28/1991 7/28/1991 3 cage_2
## 43 43 9 33 0 1 7/9/1991 4/7/1992 9 cage_2
## 44 44 3 45 1 1 12/31/1989 4/1/1990 3 cage_3
## 45 45 35 33 0 1 12/20/1989 11/18/1992 35 cage_2
## 46 46 8 28 0 1 6/22/1991 2/20/1992 8 cage_1
## 47 47 11 31 0 1 4/11/1990 3/11/1991 11 cage_2
## 48 48 56 20 1 0 5/22/1990 1/19/1995 56+ cage_1
## 49 49 2 44 0 0 11/11/1991 1/10/1992 2+ cage_3
## 50 50 3 39 1 1 1/18/1991 4/19/1991 3 cage_2
## 51 51 15 33 0 1 11/11/1989 2/10/1991 15 cage_2
## 52 52 1 31 0 1 10/1/1990 10/31/1990 1 cage_2
## 53 53 10 33 0 1 3/20/1990 1/18/1991 10 cage_2
## 54 54 1 50 1 1 7/30/1990 8/29/1990 1 cage_4
## 55 55 7 36 1 1 7/17/1989 2/14/1990 7 cage_2
## 56 56 3 30 1 1 11/10/1990 2/9/1991 3 cage_2
## 57 57 3 42 1 1 3/5/1989 6/4/1989 3 cage_3
## 58 58 2 32 1 1 3/2/1991 5/1/1991 2 cage_2
## 59 59 32 34 0 1 9/11/1989 5/11/1992 32 cage_2
## 60 60 3 38 1 1 9/12/1989 12/12/1989 3 cage_2
## 61 61 10 33 0 0 4/8/1990 2/6/1991 10+ cage_2
## 62 62 11 39 1 1 4/20/1989 3/20/1990 11 cage_2
## 63 63 3 39 1 1 1/31/1991 5/2/1991 3 cage_2
## 64 64 7 33 1 1 9/15/1989 4/15/1990 7 cage_2
## 65 65 5 34 1 1 12/7/1991 5/7/1992 5 cage_2
## 66 66 31 34 0 1 3/4/1990 10/1/1992 31 cage_2
## 67 67 5 46 1 1 4/20/1989 9/19/1989 5 cage_3
## 68 68 58 22 0 1 6/16/1989 4/15/1994 58 cage_1
## 69 69 1 44 1 1 10/1/1990 10/31/1990 1 cage_3
## 70 70 3 37 0 0 2/1/1991 5/3/1991 3+ cage_2
## 71 71 43 25 0 1 5/13/1989 12/10/1992 43 cage_1
## 72 72 1 38 0 1 8/9/1990 9/8/1990 1 cage_2
## 73 73 6 32 0 1 12/18/1991 6/17/1992 6 cage_2
## 74 74 53 34 0 1 8/23/1990 1/21/1995 53 cage_2
## 75 75 14 29 0 1 1/19/1991 3/19/1992 14 cage_1
## 76 76 4 36 1 1 8/26/1991 12/25/1991 4 cage_2
## 77 77 54 21 0 1 5/16/1991 11/13/1995 54 cage_1
## 78 78 1 26 1 1 3/20/1989 4/19/1989 1 cage_1
## 79 79 1 32 1 1 10/5/1991 11/4/1991 1 cage_2
## 80 80 8 42 0 1 5/21/1991 1/19/1992 8 cage_3
## 81 81 5 40 1 1 6/10/1991 11/9/1991 5 cage_3
## 82 82 1 37 1 1 8/31/1989 9/30/1989 1 cage_2
## 83 83 1 47 0 1 12/28/1991 1/27/1992 1 cage_3
## 84 84 2 32 1 1 9/29/1990 11/28/1990 2 cage_2
## 85 85 7 41 1 0 11/20/1991 6/19/1992 7+ cage_3
## 86 86 1 46 1 0 7/2/1989 8/1/1989 1+ cage_3
## 87 87 10 26 1 1 10/11/1991 8/10/1992 10 cage_1
## 88 88 24 30 0 0 10/11/1990 10/10/1992 24+ cage_2
## 89 89 7 32 1 1 12/5/1990 7/5/1991 7 cage_2
## 90 90 12 31 1 0 9/8/1989 9/8/1990 12+ cage_2
## 91 91 4 35 0 1 4/10/1990 8/9/1990 4 cage_2
## 92 92 57 36 0 1 12/11/1990 9/9/1995 57 cage_2
## 93 93 1 41 1 1 12/15/1990 1/14/1991 1 cage_3
## 94 94 12 36 1 0 1/13/1989 1/13/1990 12+ cage_2
## 95 95 7 35 1 1 8/22/1991 3/21/1992 7 cage_2
## 96 96 1 34 1 1 8/2/1991 9/1/1991 1 cage_2
## 97 97 5 28 0 1 5/22/1991 10/21/1991 5 cage_1
## 98 98 60 29 0 0 4/2/1990 4/1/1995 60+ cage_1
## 99 99 2 35 1 0 5/1/1991 6/30/1991 2+ cage_2
## 100 100 1 34 1 1 5/11/1989 6/10/1989 1 cage_2
cdrug <- as.factor(hmohiv$drug)
cdrug <- relevel(cdrug, ref = "1")
hmohiv$cage <- relevel(factor(hmohiv$cage), ref = "cage_3")
# Model Cox PH dengan variabel penjelas 'drug' dan 'age_cat'
coxph_drug_age <- coxph(s_obj ~ cage + cdrug, data=hmohiv,method="breslow")
summary(coxph_drug_age)
## Call:
## coxph(formula = s_obj ~ cage + cdrug, data = hmohiv, method = "breslow")
##
## n= 100, number of events= 80
##
## coef exp(coef) se(coef) z Pr(>|z|)
## cagecage_1 -1.8771 0.1530 0.4879 -3.847 0.000119 ***
## cagecage_2 -0.5991 0.5493 0.2750 -2.178 0.029385 *
## cagecage_4 1.2308 3.4239 0.6354 1.937 0.052741 .
## cdrug0 -0.8711 0.4185 0.2532 -3.440 0.000581 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## cagecage_1 0.1530 6.5345 0.05882 0.3982
## cagecage_2 0.5493 1.8204 0.32044 0.9417
## cagecage_4 3.4239 0.2921 0.98554 11.8954
## cdrug0 0.4185 2.3895 0.25479 0.6874
##
## Concordance= 0.685 (se = 0.034 )
## Likelihood ratio test= 33.11 on 4 df, p=1e-06
## Wald test = 29.43 on 4 df, p=6e-06
## Score (logrank) test = 32.93 on 4 df, p=1e-06