library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.3.2
library(carData)
## Warning: package 'carData' was built under R version 4.3.2
library(survival)
## Warning: package 'survival' was built under R version 4.3.3
library(car)
## Warning: package 'car' was built under R version 4.3.3
hmohiv<-read.table("https://stats.idre.ucla.edu/stat/r/examples/asa/hmohiv.csv", sep=",", header = TRUE)
s_obj<-Surv(hmohiv$time,hmohiv$censor)
hmohiv
## ID time age drug censor entdate enddate
## 1 1 5 46 0 1 5/15/1990 10/14/1990
## 2 2 6 35 1 0 9/19/1989 3/20/1990
## 3 3 8 30 1 1 4/21/1991 12/20/1991
## 4 4 3 30 1 1 1/3/1991 4/4/1991
## 5 5 22 36 0 1 9/18/1989 7/19/1991
## 6 6 1 32 1 0 3/18/1991 4/17/1991
## 7 7 7 36 1 1 11/11/1989 6/11/1990
## 8 8 9 31 1 1 11/25/1989 8/25/1990
## 9 9 3 48 0 1 2/11/1991 5/13/1991
## 10 10 12 47 0 1 8/11/1989 8/11/1990
## 11 11 2 28 1 0 4/11/1990 6/10/1990
## 12 12 12 34 0 1 5/11/1991 5/10/1992
## 13 13 1 44 1 1 1/17/1989 2/16/1989
## 14 14 15 32 1 1 2/16/1991 5/17/1992
## 15 15 34 36 0 1 4/9/1991 2/6/1994
## 16 16 1 36 0 1 3/9/1991 4/8/1991
## 17 17 4 54 0 1 8/3/1990 12/2/1990
## 18 18 19 35 0 0 6/10/1990 1/8/1992
## 19 19 3 44 1 0 6/12/1991 9/11/1991
## 20 20 2 38 0 1 1/7/1991 3/8/1991
## 21 21 2 40 0 0 8/29/1989 10/28/1989
## 22 22 6 34 1 1 5/29/1989 11/27/1989
## 23 23 60 25 0 0 11/16/1990 11/14/1995
## 24 24 11 32 0 1 5/9/1990 4/8/1991
## 25 25 2 42 1 0 9/10/1991 11/9/1991
## 26 26 5 47 0 1 12/26/1991 5/26/1992
## 27 27 4 30 0 0 5/29/1991 9/27/1991
## 28 28 1 47 1 1 5/1/1990 5/31/1990
## 29 29 13 41 0 1 3/24/1991 4/22/1992
## 30 30 3 40 1 1 7/18/1989 10/17/1989
## 31 31 2 43 0 1 9/16/1990 11/15/1990
## 32 32 1 41 0 1 6/22/1989 7/22/1989
## 33 33 30 30 0 1 4/27/1990 10/25/1992
## 34 34 7 37 0 1 5/16/1990 12/14/1990
## 35 35 4 42 1 1 2/19/1989 6/20/1989
## 36 36 8 31 1 1 2/17/1990 10/18/1990
## 37 37 5 39 1 1 8/6/1991 1/5/1992
## 38 38 10 32 0 1 8/10/1989 6/10/1990
## 39 39 2 51 0 1 12/27/1990 2/25/1991
## 40 40 9 36 0 1 4/26/1989 1/24/1990
## 41 41 36 43 0 1 12/4/1990 12/3/1993
## 42 42 3 39 0 1 4/28/1991 7/28/1991
## 43 43 9 33 0 1 7/9/1991 4/7/1992
## 44 44 3 45 1 1 12/31/1989 4/1/1990
## 45 45 35 33 0 1 12/20/1989 11/18/1992
## 46 46 8 28 0 1 6/22/1991 2/20/1992
## 47 47 11 31 0 1 4/11/1990 3/11/1991
## 48 48 56 20 1 0 5/22/1990 1/19/1995
## 49 49 2 44 0 0 11/11/1991 1/10/1992
## 50 50 3 39 1 1 1/18/1991 4/19/1991
## 51 51 15 33 0 1 11/11/1989 2/10/1991
## 52 52 1 31 0 1 10/1/1990 10/31/1990
## 53 53 10 33 0 1 3/20/1990 1/18/1991
## 54 54 1 50 1 1 7/30/1990 8/29/1990
## 55 55 7 36 1 1 7/17/1989 2/14/1990
## 56 56 3 30 1 1 11/10/1990 2/9/1991
## 57 57 3 42 1 1 3/5/1989 6/4/1989
## 58 58 2 32 1 1 3/2/1991 5/1/1991
## 59 59 32 34 0 1 9/11/1989 5/11/1992
## 60 60 3 38 1 1 9/12/1989 12/12/1989
## 61 61 10 33 0 0 4/8/1990 2/6/1991
## 62 62 11 39 1 1 4/20/1989 3/20/1990
## 63 63 3 39 1 1 1/31/1991 5/2/1991
## 64 64 7 33 1 1 9/15/1989 4/15/1990
## 65 65 5 34 1 1 12/7/1991 5/7/1992
## 66 66 31 34 0 1 3/4/1990 10/1/1992
## 67 67 5 46 1 1 4/20/1989 9/19/1989
## 68 68 58 22 0 1 6/16/1989 4/15/1994
## 69 69 1 44 1 1 10/1/1990 10/31/1990
## 70 70 3 37 0 0 2/1/1991 5/3/1991
## 71 71 43 25 0 1 5/13/1989 12/10/1992
## 72 72 1 38 0 1 8/9/1990 9/8/1990
## 73 73 6 32 0 1 12/18/1991 6/17/1992
## 74 74 53 34 0 1 8/23/1990 1/21/1995
## 75 75 14 29 0 1 1/19/1991 3/19/1992
## 76 76 4 36 1 1 8/26/1991 12/25/1991
## 77 77 54 21 0 1 5/16/1991 11/13/1995
## 78 78 1 26 1 1 3/20/1989 4/19/1989
## 79 79 1 32 1 1 10/5/1991 11/4/1991
## 80 80 8 42 0 1 5/21/1991 1/19/1992
## 81 81 5 40 1 1 6/10/1991 11/9/1991
## 82 82 1 37 1 1 8/31/1989 9/30/1989
## 83 83 1 47 0 1 12/28/1991 1/27/1992
## 84 84 2 32 1 1 9/29/1990 11/28/1990
## 85 85 7 41 1 0 11/20/1991 6/19/1992
## 86 86 1 46 1 0 7/2/1989 8/1/1989
## 87 87 10 26 1 1 10/11/1991 8/10/1992
## 88 88 24 30 0 0 10/11/1990 10/10/1992
## 89 89 7 32 1 1 12/5/1990 7/5/1991
## 90 90 12 31 1 0 9/8/1989 9/8/1990
## 91 91 4 35 0 1 4/10/1990 8/9/1990
## 92 92 57 36 0 1 12/11/1990 9/9/1995
## 93 93 1 41 1 1 12/15/1990 1/14/1991
## 94 94 12 36 1 0 1/13/1989 1/13/1990
## 95 95 7 35 1 1 8/22/1991 3/21/1992
## 96 96 1 34 1 1 8/2/1991 9/1/1991
## 97 97 5 28 0 1 5/22/1991 10/21/1991
## 98 98 60 29 0 0 4/2/1990 4/1/1995
## 99 99 2 35 1 0 5/1/1991 6/30/1991
## 100 100 1 34 1 1 5/11/1989 6/10/1989
- Ubah variabel ‘age’ menjadi variabel kategorik sebagai berikut,
20-29 : cage_1 30-39 : cage_2 40-49 : cage_3 50-54 : cage_4
hmohiv$cage<-recode(hmohiv$age, "20:29='cage_1'; 30:39='cage_2'; 40:49='cage_3';50:54='cage_4'", as.factor = TRUE)
hmohiv
## ID time age drug censor entdate enddate cage
## 1 1 5 46 0 1 5/15/1990 10/14/1990 cage_3
## 2 2 6 35 1 0 9/19/1989 3/20/1990 cage_2
## 3 3 8 30 1 1 4/21/1991 12/20/1991 cage_2
## 4 4 3 30 1 1 1/3/1991 4/4/1991 cage_2
## 5 5 22 36 0 1 9/18/1989 7/19/1991 cage_2
## 6 6 1 32 1 0 3/18/1991 4/17/1991 cage_2
## 7 7 7 36 1 1 11/11/1989 6/11/1990 cage_2
## 8 8 9 31 1 1 11/25/1989 8/25/1990 cage_2
## 9 9 3 48 0 1 2/11/1991 5/13/1991 cage_3
## 10 10 12 47 0 1 8/11/1989 8/11/1990 cage_3
## 11 11 2 28 1 0 4/11/1990 6/10/1990 cage_1
## 12 12 12 34 0 1 5/11/1991 5/10/1992 cage_2
## 13 13 1 44 1 1 1/17/1989 2/16/1989 cage_3
## 14 14 15 32 1 1 2/16/1991 5/17/1992 cage_2
## 15 15 34 36 0 1 4/9/1991 2/6/1994 cage_2
## 16 16 1 36 0 1 3/9/1991 4/8/1991 cage_2
## 17 17 4 54 0 1 8/3/1990 12/2/1990 cage_4
## 18 18 19 35 0 0 6/10/1990 1/8/1992 cage_2
## 19 19 3 44 1 0 6/12/1991 9/11/1991 cage_3
## 20 20 2 38 0 1 1/7/1991 3/8/1991 cage_2
## 21 21 2 40 0 0 8/29/1989 10/28/1989 cage_3
## 22 22 6 34 1 1 5/29/1989 11/27/1989 cage_2
## 23 23 60 25 0 0 11/16/1990 11/14/1995 cage_1
## 24 24 11 32 0 1 5/9/1990 4/8/1991 cage_2
## 25 25 2 42 1 0 9/10/1991 11/9/1991 cage_3
## 26 26 5 47 0 1 12/26/1991 5/26/1992 cage_3
## 27 27 4 30 0 0 5/29/1991 9/27/1991 cage_2
## 28 28 1 47 1 1 5/1/1990 5/31/1990 cage_3
## 29 29 13 41 0 1 3/24/1991 4/22/1992 cage_3
## 30 30 3 40 1 1 7/18/1989 10/17/1989 cage_3
## 31 31 2 43 0 1 9/16/1990 11/15/1990 cage_3
## 32 32 1 41 0 1 6/22/1989 7/22/1989 cage_3
## 33 33 30 30 0 1 4/27/1990 10/25/1992 cage_2
## 34 34 7 37 0 1 5/16/1990 12/14/1990 cage_2
## 35 35 4 42 1 1 2/19/1989 6/20/1989 cage_3
## 36 36 8 31 1 1 2/17/1990 10/18/1990 cage_2
## 37 37 5 39 1 1 8/6/1991 1/5/1992 cage_2
## 38 38 10 32 0 1 8/10/1989 6/10/1990 cage_2
## 39 39 2 51 0 1 12/27/1990 2/25/1991 cage_4
## 40 40 9 36 0 1 4/26/1989 1/24/1990 cage_2
## 41 41 36 43 0 1 12/4/1990 12/3/1993 cage_3
## 42 42 3 39 0 1 4/28/1991 7/28/1991 cage_2
## 43 43 9 33 0 1 7/9/1991 4/7/1992 cage_2
## 44 44 3 45 1 1 12/31/1989 4/1/1990 cage_3
## 45 45 35 33 0 1 12/20/1989 11/18/1992 cage_2
## 46 46 8 28 0 1 6/22/1991 2/20/1992 cage_1
## 47 47 11 31 0 1 4/11/1990 3/11/1991 cage_2
## 48 48 56 20 1 0 5/22/1990 1/19/1995 cage_1
## 49 49 2 44 0 0 11/11/1991 1/10/1992 cage_3
## 50 50 3 39 1 1 1/18/1991 4/19/1991 cage_2
## 51 51 15 33 0 1 11/11/1989 2/10/1991 cage_2
## 52 52 1 31 0 1 10/1/1990 10/31/1990 cage_2
## 53 53 10 33 0 1 3/20/1990 1/18/1991 cage_2
## 54 54 1 50 1 1 7/30/1990 8/29/1990 cage_4
## 55 55 7 36 1 1 7/17/1989 2/14/1990 cage_2
## 56 56 3 30 1 1 11/10/1990 2/9/1991 cage_2
## 57 57 3 42 1 1 3/5/1989 6/4/1989 cage_3
## 58 58 2 32 1 1 3/2/1991 5/1/1991 cage_2
## 59 59 32 34 0 1 9/11/1989 5/11/1992 cage_2
## 60 60 3 38 1 1 9/12/1989 12/12/1989 cage_2
## 61 61 10 33 0 0 4/8/1990 2/6/1991 cage_2
## 62 62 11 39 1 1 4/20/1989 3/20/1990 cage_2
## 63 63 3 39 1 1 1/31/1991 5/2/1991 cage_2
## 64 64 7 33 1 1 9/15/1989 4/15/1990 cage_2
## 65 65 5 34 1 1 12/7/1991 5/7/1992 cage_2
## 66 66 31 34 0 1 3/4/1990 10/1/1992 cage_2
## 67 67 5 46 1 1 4/20/1989 9/19/1989 cage_3
## 68 68 58 22 0 1 6/16/1989 4/15/1994 cage_1
## 69 69 1 44 1 1 10/1/1990 10/31/1990 cage_3
## 70 70 3 37 0 0 2/1/1991 5/3/1991 cage_2
## 71 71 43 25 0 1 5/13/1989 12/10/1992 cage_1
## 72 72 1 38 0 1 8/9/1990 9/8/1990 cage_2
## 73 73 6 32 0 1 12/18/1991 6/17/1992 cage_2
## 74 74 53 34 0 1 8/23/1990 1/21/1995 cage_2
## 75 75 14 29 0 1 1/19/1991 3/19/1992 cage_1
## 76 76 4 36 1 1 8/26/1991 12/25/1991 cage_2
## 77 77 54 21 0 1 5/16/1991 11/13/1995 cage_1
## 78 78 1 26 1 1 3/20/1989 4/19/1989 cage_1
## 79 79 1 32 1 1 10/5/1991 11/4/1991 cage_2
## 80 80 8 42 0 1 5/21/1991 1/19/1992 cage_3
## 81 81 5 40 1 1 6/10/1991 11/9/1991 cage_3
## 82 82 1 37 1 1 8/31/1989 9/30/1989 cage_2
## 83 83 1 47 0 1 12/28/1991 1/27/1992 cage_3
## 84 84 2 32 1 1 9/29/1990 11/28/1990 cage_2
## 85 85 7 41 1 0 11/20/1991 6/19/1992 cage_3
## 86 86 1 46 1 0 7/2/1989 8/1/1989 cage_3
## 87 87 10 26 1 1 10/11/1991 8/10/1992 cage_1
## 88 88 24 30 0 0 10/11/1990 10/10/1992 cage_2
## 89 89 7 32 1 1 12/5/1990 7/5/1991 cage_2
## 90 90 12 31 1 0 9/8/1989 9/8/1990 cage_2
## 91 91 4 35 0 1 4/10/1990 8/9/1990 cage_2
## 92 92 57 36 0 1 12/11/1990 9/9/1995 cage_2
## 93 93 1 41 1 1 12/15/1990 1/14/1991 cage_3
## 94 94 12 36 1 0 1/13/1989 1/13/1990 cage_2
## 95 95 7 35 1 1 8/22/1991 3/21/1992 cage_2
## 96 96 1 34 1 1 8/2/1991 9/1/1991 cage_2
## 97 97 5 28 0 1 5/22/1991 10/21/1991 cage_1
## 98 98 60 29 0 0 4/2/1990 4/1/1995 cage_1
## 99 99 2 35 1 0 5/1/1991 6/30/1991 cage_2
## 100 100 1 34 1 1 5/11/1989 6/10/1989 cage_2
cage.coxph <- coxph( s_obj ~ hmohiv$cage, data = hmohiv, method="breslow")
summary(cage.coxph)
## Call:
## coxph(formula = s_obj ~ hmohiv$cage, data = hmohiv, method = "breslow")
##
## n= 100, number of events= 80
##
## coef exp(coef) se(coef) z Pr(>|z|)
## hmohiv$cagecage_2 1.2552 3.5086 0.4315 2.909 0.00362 **
## hmohiv$cagecage_3 1.7974 6.0340 0.4762 3.774 0.00016 ***
## hmohiv$cagecage_4 2.7635 15.8552 0.7350 3.760 0.00017 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## hmohiv$cagecage_2 3.509 0.28501 1.506 8.174
## hmohiv$cagecage_3 6.034 0.16573 2.373 15.345
## hmohiv$cagecage_4 15.855 0.06307 3.755 66.954
##
## Concordance= 0.63 (se = 0.033 )
## Likelihood ratio test= 21.29 on 3 df, p=9e-05
## Wald test = 19.28 on 3 df, p=2e-04
## Score (logrank) test = 21.96 on 3 df, p=7e-05
- Buat model Cox PH dengan variabel penjelas drug (ref:drug=1) dan
cage (ref: cage_3)
cdrug <- as.factor(hmohiv$drug)
cdrug <- relevel(cdrug, ref = "1")
hmohiv$cage <- relevel(hmohiv$cage, ref = "cage_3")
cox_model <- coxph(s_obj ~ cage + drug , data = hmohiv, method = "breslow")
summary(cox_model)
## Call:
## coxph(formula = s_obj ~ cage + drug, 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 .
## drug 0.8711 2.3895 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
## drug 2.3895 0.4185 1.45475 3.9247
##
## 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
- Interpretasikan hasilnya!