statistic-celiac-doc
Zia Farooq
2022-07-04
library(survival)
library(survminer)## Loading required package: ggplot2## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomalibrary(simstudy) #the package for non-PH models
library(ggplot2)
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, unionlibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ tibble 3.1.7 ✔ purrr 0.3.4
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(viridis)## Loading required package: viridisLitelibrary(ggrepel)
library(ggthemes)
library(lubridate) # to handle the dates data##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(anytime) #handles the dates data in any format
library(fmsb) # for 95% CI of incidence rateUpload the data files
celiac <- read.csv("~/Desktop/Celiac-main/celiac.txt")
nonceliac <- read.csv("~/Desktop/Celiac-main/nonceliac.txt")
ICD <- read.csv("~/Desktop/Celiac-main/ICD.txt")
# peek into the data
head(celiac)## id.number Start.of.flwp Gender Age Birthday Date.of.Death Date.of.Migration
## 1 424 20081106 2 72 1936-05-15 NA NA
## 2 881 20070613 1 21 1986-02-15 NA NA
## 3 1739 20040609 2 12 1991-12-15 NA NA
## 4 2462 20000306 2 19 1980-09-15 NA NA
## 5 2652 20040712 2 46 1957-10-15 NA NA
## 6 2785 20031001 1 17 1986-01-15 NA NAhead(nonceliac)## id.number Celiac.number Date.of.Death Date.of.Migration
## 1 7659518 424 NA NA
## 2 1038521 424 NA NA
## 3 9265189 881 NA NA
## 4 10658482 881 NA NA
## 5 8460019 1739 NA NA
## 6 10980993 1739 NA NAhead(ICD)## id.number ICD.codes Date.of.ICD.codes
## 1 2855591 C402 2003-03-04
## 2 2855591 C402 2002-12-25
## 3 5552533 C787 C809 M069 2010-05-25
## 4 9515567 D489 C921 2001-09-20
## 5 2112819 Z489 C153 2002-09-11
## 6 9075074 C900 2005-09-23Data cleaning step
Dates columns need to be in one format
celiac <- celiac %>%
mutate(across(c(5:7), as.character))%>%
mutate(across(c(2,5,6,7), anydate))
nonceliac <- nonceliac %>%
mutate(across(c(3,4), as.character))%>%
mutate(across(c(3,4), anydate))
ICD <-ICD%>%
mutate(Date.of.ICD.codes= as.Date(Date.of.ICD.codes))
#check the data sets
head(celiac)## id.number Start.of.flwp Gender Age Birthday Date.of.Death Date.of.Migration
## 1 424 2008-11-06 2 72 1936-05-15 <NA> <NA>
## 2 881 2007-06-13 1 21 1986-02-15 <NA> <NA>
## 3 1739 2004-06-09 2 12 1991-12-15 <NA> <NA>
## 4 2462 2000-03-06 2 19 1980-09-15 <NA> <NA>
## 5 2652 2004-07-12 2 46 1957-10-15 <NA> <NA>
## 6 2785 2003-10-01 1 17 1986-01-15 <NA> <NA>head(nonceliac)## id.number Celiac.number Date.of.Death Date.of.Migration
## 1 7659518 424 <NA> <NA>
## 2 1038521 424 <NA> <NA>
## 3 9265189 881 <NA> <NA>
## 4 10658482 881 <NA> <NA>
## 5 8460019 1739 <NA> <NA>
## 6 10980993 1739 <NA> <NA>head(ICD)## id.number ICD.codes Date.of.ICD.codes
## 1 2855591 C402 2003-03-04
## 2 2855591 C402 2002-12-25
## 3 5552533 C787 C809 M069 2010-05-25
## 4 9515567 D489 C921 2001-09-20
## 5 2112819 Z489 C153 2002-09-11
## 6 9075074 C900 2005-09-23Extract CRC events data from the ICD file
# Find out the CRC events from ICD file both for celiac and nonceliac
ICD_CRC_events <- ICD %>%
filter(grepl('C18', ICD.codes))%>%
rbind(
ICD %>% filter(grepl('C19', ICD.codes)) )%>%
rbind(
ICD %>% filter(grepl('C20', ICD.codes)) )%>%
mutate(Date.of.ICD.codes= as.Date(Date.of.ICD.codes))
head(ICD_CRC_events)## id.number ICD.codes Date.of.ICD.codes
## 1 10447614 C187 C787 C771 C722 I109 E119 2009-04-20
## 2 3912733 C187 2005-01-24
## 3 2040543 C189 2005-01-28
## 4 4780187 C186 I109 E785 K900A 2010-08-25
## 5 344621 C186 I489B I209 N390 2008-06-07
## 6 7462794 Z514 C183 2010-01-06Construct the data frame for celiac patients based on the instructions
#Find first CRC events both for celiac and nonceliac groups
celiac_CRC <- celiac%>%
filter(id.number %in% unique(ICD_CRC_events$id.number) ) %>%
merge(ICD_CRC_events%>%filter(id.number %in% unique(celiac$id.number))%>%
group_by(id.number)%>% summarise(Date.of.CRC.onset = min(Date.of.ICD.codes) ),
by.x = c("id.number"),
by.y =c("id.number"),
all.x = TRUE)%>%
filter(Start.of.flwp <= Date.of.CRC.onset)%>% # exclude the CRC cases before the follow-up
mutate(Date.of.flwp.end = do.call(pmin, c(select(., c(6:8)),
na.rm = TRUE)) )%>% # Case-wise End date of follow-up conditioned on Death/Migration/CRC case
#filter(Start.of.flwp>Date.of.flwp.end)%>%
mutate(person.time.years = as.numeric(Date.of.CRC.onset-Start.of.flwp)/ 365.25)%>% # calculate person-years time
mutate(CRC.status= "Yes") #mutate the variable for the CRC events
celiac_nonCRC <- celiac %>%filter(id.number %in% setdiff(unique(celiac$id.number),unique(ICD_CRC_events$id.number) ))%>%
mutate(Date.of.CRC.onset= as.Date(NA))%>%
mutate(Date.of.flwp.end = as.Date("2019-12-31") )%>%
#filter(Start.of.flwp <= Date.of.CRC.onset)%>%
mutate(Date.of.flwp.end = do.call(pmin, c(select(., c(6:9)),
na.rm = TRUE)) )%>%
#filter(Start.of.flwp>Date.of.flwp.end)%>%
mutate(person.time.years = as.numeric(Date.of.flwp.end-Start.of.flwp)/ 365.25)%>%
mutate(CRC.status= "No") #mutate the variable for the CRC events
celiac_final_df <- rbind(celiac_CRC, celiac_nonCRC)%>%
filter(Start.of.flwp<Date.of.flwp.end)%>% # filter out the events with negative follow-up time
mutate(CRC.status= as.factor(CRC.status))%>%
mutate(exposure.status= "1") #mutate the variable for the exposure status (celiac~1, nonceliac~0)
dim(celiac_CRC)## [1] 113 11dim(celiac_nonCRC)## [1] 12199 11dim(celiac_final_df)## [1] 12310 12Calculate the required estimates (events,person years time/1000 and Incidence rate) for Celiac patients
total_cases_celiac= length(celiac_CRC$id.number)
person_time_celiac= sum(celiac_final_df$person.time.years)
Incidence_rate_celiac = total_cases_celiac/ person_time_celiac *1000
total_cases_celiac## [1] 113person_time_celiac## [1] 137267.8## 95% CI for incidence rate using fmsb package
IRCI(total_cases_celiac*1000, person_time_celiac, conf.level=0.95) ## $IR
## [1] 0.8232082
##
## $IRL
## [1] 0.8184085
##
## $IRU
## [1] 0.828008# OR using the following lines
SE_celiac = 1.96*sqrt((total_cases_celiac*1000)/person_time_celiac^2 )
CIL_95_celiac_IR <- Incidence_rate_celiac-SE_celiac
CIU_95_celiac_IR <- Incidence_rate_celiac+SE_celiacConstruct the data frame for nonceliac (nonexposed) group based on the instructions
nonceliac_CRC <-nonceliac%>%merge(celiac[,c(1:5)], by.x = "Celiac.number" , by.y = "id.number" ,all.x = TRUE )%>%
filter(id.number %in% unique(ICD_CRC_events$id.number) ) %>%
merge( ICD_CRC_events%>%filter(id.number %in% unique(nonceliac$id.number)) %>%
group_by(id.number)%>% summarise(Date.of.CRC.onset = min(Date.of.ICD.codes)),
by.x = c("id.number"),
by.y =c("id.number"),
all.y = TRUE
)%>%
select(c(1,5:8,3,4,9))%>%
filter(Start.of.flwp <= Date.of.CRC.onset)%>%
mutate(Date.of.flwp.end = do.call(pmin, c(select(., c(6:8)),
na.rm = TRUE)) )%>%
filter(Start.of.flwp<Date.of.flwp.end)%>%
mutate(person.time.years = as.numeric(Date.of.CRC.onset-Start.of.flwp)/ 365.25)%>%
mutate(CRC.status= "Yes")
nonceliac_nonCRC <- nonceliac%>%
#match the nonceliac patients with celiac with celiac number
merge(celiac[,c(1:5)], by.x = "Celiac.number" , by.y = "id.number" ,all.x = TRUE )%>%
#filter the nonceliacs found in the ICD data
filter(id.number %in% setdiff(unique(nonceliac$id.number),unique(ICD_CRC_events$id.number) ))%>%
#mutate 'Date.of.CRC.onset= NA' variable
mutate(Date.of.CRC.onset= as.Date(NA))%>%
#mutate Followup end date variable with fixed date(2019-12-31)
mutate(Date.of.flwp.end = as.Date("2019-12-31") )%>%
select(c(2,5:8,3:4,9,10))%>%
#
mutate(Date.of.flwp.end = do.call(pmin, c(select(., c(6:9)),
na.rm = TRUE)) )%>%
filter(Start.of.flwp<Date.of.flwp.end)%>%
mutate(person.time.years = as.numeric(Date.of.flwp.end-Start.of.flwp)/ 365.25)%>%
mutate(CRC.status= "No")
nonceliac_final_df <- rbind(nonceliac_CRC, nonceliac_nonCRC)%>%
filter(Start.of.flwp<Date.of.flwp.end)%>%
mutate(CRC.status= as.factor(CRC.status))%>%
mutate(exposure.status= "0")Calculate the required estimates (events,person years time /1000 and Incidence rate) for nonceliac persons
total_cases_nonceliac= length(nonceliac_CRC$id.number)
person_time_nonceliac= sum(nonceliac_final_df$person.time.years)
Incidence_rate_nonceliac = (total_cases_nonceliac/ person_time_nonceliac *1000)
total_cases_nonceliac## [1] 194person_time_nonceliac## [1] 272313.7## 95% CI for incidence rate using fmsb package
IRCI(total_cases_nonceliac*1000, person_time_nonceliac, conf.level=0.95) ## $IR
## [1] 0.7124136
##
## $IRL
## [1] 0.7092435
##
## $IRU
## [1] 0.7155838# OR
SE_nonceliac = 1.96*sqrt((total_cases_nonceliac*1000)/person_time_nonceliac^2 )
CIL_95_nonceliac_IR <- Incidence_rate_nonceliac-SE_nonceliac
CIL_95_nonceliac_IR <- Incidence_rate_nonceliac+SE_nonceliacConstruct the final data frame to estimate the hazard ratios (HR) using Cox proportion hazard model
df_celiac_nonceliac <- rbind(celiac_final_df,nonceliac_final_df)%>%
mutate(flwp.time =as.numeric((Date.of.flwp.end-Start.of.flwp)) )%>%
mutate(exposure.status=as.factor(exposure.status))%>%
mutate(Gender= as.factor(Gender))%>%
filter(Start.of.flwp<Date.of.flwp.end)%>%
select(c(1:10,12,11,13))
head(df_celiac_nonceliac)## id.number Start.of.flwp Gender Age Birthday Date.of.Death Date.of.Migration
## 1 18671 2002-02-01 1 77 1924-11-15 2003-08-27 <NA>
## 2 41943 2000-03-22 2 74 1925-08-15 2004-04-15 <NA>
## 3 170391 2000-05-03 1 37 1963-01-15 2017-07-09 <NA>
## 4 183428 2005-06-01 2 34 1970-11-15 <NA> <NA>
## 5 199935 2005-03-23 1 68 1937-02-15 2006-01-01 <NA>
## 6 224099 2010-06-09 1 76 1933-08-15 2015-02-03 <NA>
## Date.of.CRC.onset Date.of.flwp.end person.time.years exposure.status
## 1 2002-04-09 2002-04-09 0.1834360 1
## 2 2003-04-21 2003-04-21 3.0800821 1
## 3 2015-08-31 2015-08-31 15.3264887 1
## 4 2014-04-29 2014-04-29 8.9089665 1
## 5 2005-08-21 2005-08-21 0.4134155 1
## 6 2013-01-21 2013-01-21 2.6201232 1
## CRC.status flwp.time
## 1 Yes 67
## 2 Yes 1125
## 3 Yes 5598
## 4 Yes 3254
## 5 Yes 151
## 6 Yes 957Apply Cox hazard model for unadjusted
HR_CRC_Unadjusted <- coxph(Surv( as.numeric(Start.of.flwp),
as.numeric(Date.of.flwp.end), as.numeric(CRC.status))
~exposure.status, data=df_celiac_nonceliac)
summary(HR_CRC_Unadjusted)## Call:
## coxph(formula = Surv(as.numeric(Start.of.flwp), as.numeric(Date.of.flwp.end),
## as.numeric(CRC.status)) ~ exposure.status, data = df_celiac_nonceliac)
##
## n= 36668, number of events= 307
##
## coef exp(coef) se(coef) z Pr(>|z|)
## exposure.status1 0.1443 1.1552 0.1183 1.219 0.223
##
## exp(coef) exp(-coef) lower .95 upper .95
## exposure.status1 1.155 0.8657 0.916 1.457
##
## Concordance= 0.497 (se = 0.014 )
## Likelihood ratio test= 1.47 on 1 df, p=0.2
## Wald test = 1.49 on 1 df, p=0.2
## Score (logrank) test = 1.49 on 1 df, p=0.2#Test PH assumption
test.ph <- cox.zph(HR_CRC_Unadjusted)
#test.ph
ggcoxzph(test.ph)## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
Apply the Cox model adjusted for Age and Gender variables
HR_CRC_AgeGender_adjusted <- coxph(Surv(flwp.time,as.numeric(CRC.status) )
~exposure.status+Age+Gender, data=df_celiac_nonceliac)
summary(HR_CRC_AgeGender_adjusted)## Call:
## coxph(formula = Surv(flwp.time, as.numeric(CRC.status)) ~ exposure.status +
## Age + Gender, data = df_celiac_nonceliac)
##
## n= 36668, number of events= 307
##
## coef exp(coef) se(coef) z Pr(>|z|)
## exposure.status1 0.135881 1.145546 0.118374 1.148 0.251
## Age 0.072753 1.075465 0.003578 20.332 <2e-16 ***
## Gender2 -0.154423 0.856910 0.114902 -1.344 0.179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## exposure.status1 1.1455 0.8729 0.9083 1.445
## Age 1.0755 0.9298 1.0679 1.083
## Gender2 0.8569 1.1670 0.6841 1.073
##
## Concordance= 0.888 (se = 0.006 )
## Likelihood ratio test= 668.5 on 3 df, p=<2e-16
## Wald test = 420.9 on 3 df, p=<2e-16
## Score (logrank) test = 665.8 on 3 df, p=<2e-16PH_assump_CRC_AgeGender_Adjusted <- cox.zph(HR_CRC_AgeGender_adjusted)
PH_assump_CRC_AgeGender_Adjusted## chisq df p
## exposure.status 21.949 1 2.8e-06
## Age 2.814 1 0.093
## Gender 0.489 1 0.484
## GLOBAL 24.748 3 1.7e-05ggcoxzph(PH_assump_CRC_AgeGender_Adjusted)
Stratify the followup duration time
#Divide the followup duration into 3 intervals [0,200), [200, 1650), [1650, 7302) #Use surveSplit function from simstudy package
df_celiac_nonceliac_stratfied <- survSplit(Surv(flwp.time, as.numeric(CRC.status)) ~.,
data= df_celiac_nonceliac, cut=c(200,1650),
episode= "tgroup")
head(df_celiac_nonceliac_stratfied)## id.number Start.of.flwp Gender Age Birthday Date.of.Death Date.of.Migration
## 1 18671 2002-02-01 1 77 1924-11-15 2003-08-27 <NA>
## 2 41943 2000-03-22 2 74 1925-08-15 2004-04-15 <NA>
## 3 41943 2000-03-22 2 74 1925-08-15 2004-04-15 <NA>
## 4 170391 2000-05-03 1 37 1963-01-15 2017-07-09 <NA>
## 5 170391 2000-05-03 1 37 1963-01-15 2017-07-09 <NA>
## 6 170391 2000-05-03 1 37 1963-01-15 2017-07-09 <NA>
## Date.of.CRC.onset Date.of.flwp.end person.time.years exposure.status tstart
## 1 2002-04-09 2002-04-09 0.183436 1 0
## 2 2003-04-21 2003-04-21 3.080082 1 0
## 3 2003-04-21 2003-04-21 3.080082 1 200
## 4 2015-08-31 2015-08-31 15.326489 1 0
## 5 2015-08-31 2015-08-31 15.326489 1 200
## 6 2015-08-31 2015-08-31 15.326489 1 1650
## flwp.time event tgroup
## 1 67 1 1
## 2 200 0 1
## 3 1125 1 2
## 4 200 0 1
## 5 1650 0 2
## 6 5598 1 3Apply again the Cox PH model on stratified followup duration- Unadjusted
HR_CRC_Unadjusted_stratfied <- survival::coxph(Surv(tstart,flwp.time, event)~
exposure.status*strata(tgroup),
data=df_celiac_nonceliac_stratfied)
summary(HR_CRC_Unadjusted_stratfied)## Call:
## survival::coxph(formula = Surv(tstart, flwp.time, event) ~ exposure.status *
## strata(tgroup), data = df_celiac_nonceliac_stratfied)
##
## n= 106558, number of events= 307
##
## coef exp(coef) se(coef) z
## exposure.status1 2.69983 14.87716 0.53229 5.072
## exposure.status1:strata(tgroup)tgroup=2 -2.50733 0.08149 0.56974 -4.401
## exposure.status1:strata(tgroup)tgroup=3 -3.06894 0.04647 0.56035 -5.477
## Pr(>|z|)
## exposure.status1 3.93e-07 ***
## exposure.status1:strata(tgroup)tgroup=2 1.08e-05 ***
## exposure.status1:strata(tgroup)tgroup=3 4.33e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## exposure.status1 14.87716 0.06722 5.24123
## exposure.status1:strata(tgroup)tgroup=2 0.08149 12.27207 0.02668
## exposure.status1:strata(tgroup)tgroup=3 0.04647 21.51909 0.01550
## upper .95
## exposure.status1 42.2286
## exposure.status1:strata(tgroup)tgroup=2 0.2489
## exposure.status1:strata(tgroup)tgroup=3 0.1394
##
## Concordance= 0.563 (se = 0.015 )
## Likelihood ratio test= 49.8 on 3 df, p=9e-11
## Wald test = 31.07 on 3 df, p=8e-07
## Score (logrank) test = 51.08 on 3 df, p=5e-11cox.zph(HR_CRC_Unadjusted_stratfied)## chisq df p
## exposure.status 0.616 1 0.43
## exposure.status:strata(tgroup) 0.570 2 0.75
## GLOBAL 5.712 3 0.13#ggcoxzph(HR_CRC_Unadjusted_stratfied)Apply again the Cox PH model on stratified followup duration- Adjusted for Age and Gender variables
HR_CRC_AgeGender_adjusted_stratified <- survival::coxph(Surv(tstart,flwp.time, event)~
exposure.status*strata(tgroup)+Age+Gender,
data=df_celiac_nonceliac_stratfied)
summary(HR_CRC_AgeGender_adjusted_stratified)## Call:
## survival::coxph(formula = Surv(tstart, flwp.time, event) ~ exposure.status *
## strata(tgroup) + Age + Gender, data = df_celiac_nonceliac_stratfied)
##
## n= 106558, number of events= 307
##
## coef exp(coef) se(coef) z
## exposure.status1 2.664553 14.361534 0.532301 5.006
## Age 0.072778 1.075492 0.003593 20.254
## Gender2 -0.153608 0.857608 0.114911 -1.337
## exposure.status1:strata(tgroup)tgroup=2 -2.480087 0.083736 0.569740 -4.353
## exposure.status1:strata(tgroup)tgroup=3 -3.037078 0.047975 0.560359 -5.420
## Pr(>|z|)
## exposure.status1 5.57e-07 ***
## Age < 2e-16 ***
## Gender2 0.181
## exposure.status1:strata(tgroup)tgroup=2 1.34e-05 ***
## exposure.status1:strata(tgroup)tgroup=3 5.96e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## exposure.status1 14.36153 0.06963 5.05947
## Age 1.07549 0.92981 1.06794
## Gender2 0.85761 1.16603 0.68466
## exposure.status1:strata(tgroup)tgroup=2 0.08374 11.94230 0.02741
## exposure.status1:strata(tgroup)tgroup=3 0.04797 20.84425 0.01600
## upper .95
## exposure.status1 40.7658
## Age 1.0831
## Gender2 1.0742
## exposure.status1:strata(tgroup)tgroup=2 0.2558
## exposure.status1:strata(tgroup)tgroup=3 0.1439
##
## Concordance= 0.896 (se = 0.006 )
## Likelihood ratio test= 715.7 on 5 df, p=<2e-16
## Wald test = 445.5 on 5 df, p=<2e-16
## Score (logrank) test = 714.5 on 5 df, p=<2e-16cox.zph(HR_CRC_AgeGender_adjusted_stratified)## chisq df p
## exposure.status 0.623 1 0.43
## Age 2.464 1 0.12
## Gender 0.453 1 0.50
## exposure.status:strata(tgroup) 0.570 2 0.75
## GLOBAL 8.436 5 0.13#ggcoxzph(HR_CRC_AgeGender_adjusted_stratified)