# Hodgkin's Disease
# Dr. Thirachit (cthirachit@yahoo.com)
# Dr. Shevachut (karshevachut@hotmail.com)
options(scipen=99999)
options("width"=200)
setwd("D:/Consults/Thirachit")
library(readxl)
library(epicalc)
## Warning: package 'epicalc' was built under R version 3.5.2
## Loading required package: foreign
## Loading required package: survival
## Loading required package: MASS
## Loading required package: nnet
zap()
noyes <- c("No","Yes")
source("summ.R")
data <- as.data.frame(read_excel("Hodgkins.xlsx"))
use(data)
# Part I: Patient characteristics ####
# case
any(duplicated(case))
## [1] FALSE
# 13 digit id
any(duplicated(id))
## [1] TRUE
sum(duplicated(id))
## [1] 11
sum(is.na(id))
## [1] 12
# Hospital number
any(duplicated(hn))
## [1] FALSE
sex <- factor(sex, labels=c("Male","Female"))
label.var(sex, "Sex")
tab1(sex, graph=FALSE)
## sex : Sex
## Frequency Percent Cum. percent
## Male 57 79.2 79.2
## Female 15 20.8 100.0
## Total 72 100.0 100.0
# Place of birth
songkhla <- factor(ifelse(prov==1, "Songkhla","Other"))
label.var(songkhla, "Born in Songkhla")
birthpv <- factor(prov, labels=c("Songkhla","Phattalung","Trang","Yala","Pattani","Narathiwat","Satun","Nakon Sri Thammarat","Surat Thani","Ranong","Phuket","Krabi","Pha-nga","Chumphon","Other"), levels=1:15)
label.var(birthpv, "Place of birth")
tab1(birthpv, graph=FALSE)
## birthpv : Place of birth
## Frequency Percent Cum. percent
## Songkhla 13 18.1 18.1
## Phattalung 2 2.8 20.8
## Trang 6 8.3 29.2
## Yala 6 8.3 37.5
## Pattani 9 12.5 50.0
## Narathiwat 8 11.1 61.1
## Satun 0 0.0 61.1
## Nakon Sri Thammarat 14 19.4 80.6
## Surat Thani 5 6.9 87.5
## Ranong 0 0.0 87.5
## Phuket 4 5.6 93.1
## Krabi 4 5.6 98.6
## Pha-nga 1 1.4 100.0
## Chumphon 0 0.0 100.0
## Other 0 0.0 100.0
## Total 72 100.0 100.0
# Date of birth
bd <- as.Date(bd)
label.var(bd, "Date of birth")
summ(bd)
## obs. mean median s.d. min. max.
## 72 2538-10-16 2541-02-04 <NA> 2515-03-29 2554-10-14
# Date of first visit
n1st <- as.Date(n1st)
label.var(n1st, "Date of first visit")
summ(n1st)
## obs. mean median s.d. min. max.
## 72 2547-08-28 2550-05-18 <NA> 2528-03-27 2560-08-17
# Date of diagnosis
datedx <- as.Date(datedx)
label.var(datedx, "Date of diagnosis")
summ(datedx)
## obs. mean median s.d. min. max.
## 72 2547-09-10 2550-05-18 <NA> 2528-03-27 2560-08-21
age2 <- as.integer(difftime(datedx, bd))/365
label.var(age2, "Age at diagnosis")
summ(age2) # Median age = 8 years (range: 3.5 - 15)
## obs. mean median s.d. min. max.
## 72 8.906 8.012 3.277 3.548 15.11
summ(age)
## obs. mean median s.d. min. max.
## 72 106.175 96 39.442 36 177.9
summ(age/12-age2) # Looks ok now
## obs. mean median s.d. min. max.
## 72 -0.058 -0.07 0.325 -0.712 0.655
# Age group
agegp <- cut(age2, breaks=c(0,5.99,10.99,99), labels=c("0-5","6-10","10-15"))
tab1(agegp, graph=FALSE)
## agegp :
## Frequency Percent Cum. percent
## 0-5 16 22.2 22.2
## 6-10 34 47.2 69.4
## 10-15 22 30.6 100.0
## Total 72 100.0 100.0
# Median age
agem <- factor(ifelse(age<median(age), "Young","Old"))
label.var(agem, "Median age")
tab1(agem, graph=FALSE)
## agem : Median age
## Frequency Percent Cum. percent
## Old 38 52.8 52.8
## Young 34 47.2 100.0
## Total 72 100.0 100.0
# Body weight at first visit (kg)
summ(wt)
## obs. mean median s.d. min. max.
## 72 26.664 22.65 13.097 12 88
plot(age/12, wt, xlab="Age (years)")
plot(ht, wt, xlab="Height (cm)", ylab="Weight (kg)")
with(.data[wt>80,],
text(ht, wt, labels=case, pos=4)
)
data.frame(case, hn, sex, age, ht, wt)[wt>80,] # Confirmed
## case hn sex age ht wt
## 61 20 962943 Female 144 158 88
# height at first visit (cm)
summ(ht)
## obs. mean median s.d. min. max.
## 72 125.642 121.25 18.009 95.5 166.5
# Monthly family income
summ(income, graph=FALSE)
## obs. mean median s.d. min. max.
## 8 16500 17500 10460.811 5000 37000
# Child sequence in family
tab1(child, graph=FALSE)
## child :
## Frequency Percent Cum. percent
## 1 7 9.7 9.7
## 2 27 37.5 47.2
## 3 19 26.4 73.6
## 4 9 12.5 86.1
## 5 3 4.2 90.3
## 6 2 2.8 93.1
## 7 1 1.4 94.4
## 8 2 2.8 97.2
## 9 1 1.4 98.6
## 12 1 1.4 100.0
## Total 72 100.0 100.0
# Insurance status
tab1(`if`, graph=FALSE)
## if :
## Frequency %(NA+) %(NA-)
## immigrant worker 1 1.4 100
## <NA> 71 98.6 0
## Total 72 100.0 100
right <- factor(right, labels=c("UC","Govt","Immigrant"))
label.var(right, "Type of insurance")
tab1(right, graph=FALSE)
## right : Type of insurance
## Frequency Percent Cum. percent
## UC 67 93.1 93.1
## Govt 4 5.6 98.6
## Immigrant 1 1.4 100.0
## Total 72 100.0 100.0
# Underlying disease
ud <- factor(ud, labels=noyes, levels=c(2,1))
label.var(ud, "Underlying disease")
tab1(ud, graph=FALSE)
## ud : Underlying disease
## Frequency Percent Cum. percent
## No 68 94.4 94.4
## Yes 4 5.6 100.0
## Total 72 100.0 100.0
tab1(if1, graph=FALSE)
## if1 :
## Frequency %(NA+) %(NA-)
## allergic rhinitis 1 1.4 25
## asthma 2 2.8 50
## Thalassemia disease (Hb H disease) 1 1.4 25
## <NA> 68 94.4 0
## Total 72 100.0 100
# Part II: Clinical presentation ####
# B Symptoms
fever <- factor(fever, labels=noyes, levels=c(2,1))
label.var(fever, "Fever")
tab1(fever, graph=FALSE)
## fever : Fever
## Frequency Percent Cum. percent
## No 47 65.3 65.3
## Yes 25 34.7 100.0
## Total 72 100.0 100.0
sweat <- factor(sweat, labels=noyes, levels=c(2,1))
label.var(sweat, "Night sweating")
tab1(sweat, graph=FALSE)
## sweat : Night sweating
## Frequency Percent Cum. percent
## No 68 94.4 94.4
## Yes 4 5.6 100.0
## Total 72 100.0 100.0
wtloss <- factor(wtloss, labels=noyes, levels=c(2,1))
label.var(wtloss, "Weight loss")
tab1(wtloss, graph=FALSE)
## wtloss : Weight loss
## Frequency Percent Cum. percent
## No 59 81.9 81.9
## Yes 13 18.1 100.0
## Total 72 100.0 100.0
# Any B symptoms (fever, sweat, weight loss)
bsympt <- factor(ifelse(fever=="Yes" | sweat=="Yes" | wtloss=="Yes", "Yes", "No"))
label.var(bsympt, "B symptoms")
tab1(bsympt, graph=FALSE)
## bsympt : B symptoms
## Frequency Percent Cum. percent
## No 40 55.6 55.6
## Yes 32 44.4 100.0
## Total 72 100.0 100.0
# Lymphadenopathy (1=Yes, 2=No)
label.var(ln1, "Mediastinal")
label.var(ln2, "Intra-abdominal")
label.var(ln3, "Intra-thoracic")
label.var(ln4, "Axillary")
label.var(ln5, "Inguinal")
label.var(ln6, "Pelvic")
label.var(ln7, "Cervical")
label.var(ln8, "Epitrochlear")
label.var(ln9, "Other sites")
tableStack(c(ln1:ln8, ln9))
## 1 2 count mean sd description
## ln1 19 53 72 1.7 0.4 Mediastinal
## ln2 32 40 72 1.6 0.5 Intra-abdominal
## ln3 21 51 72 1.7 0.5 Intra-thoracic
## ln4 21 51 72 1.7 0.5 Axillary
## ln5 19 53 72 1.7 0.4 Inguinal
## ln6 0 72 72 2 0 Pelvic
## ln7 60 12 72 1.2 0.4 Cervical
## ln8 1 71 72 2 0.1 Epitrochlear
## ln9 17 55 72 1.8 0.4 Other sites
## Total score 72 15.4 1.3
## Average score 72 1.7 0.1
ln8sz[!is.na(ln8sz)]
## [1] "epitrochlear node"
recode.is.na(if2, "_None_")
label.var(if2, "Other lymphadenopathy")
tab1(if2, graph=FALSE, sort="decreasing")
## if2 : Other lymphadenopathy
## Frequency Percent Cum. percent
## _None_ 55 76.4 76.4
## abdominal mass 7 9.7 86.1
## mediastinal mass 5 6.9 93.1
## splenic mass 1 1.4 94.4
## pericardial effusion, pleural effusion 1 1.4 95.8
## paravertebral mass 1 1.4 97.2
## mediastinal mass, hepatic mass 1 1.4 98.6
## chest wall mass 1 1.4 100.0
## Total 72 100.0 100.0
# Number of lymph node sites
recode(c(ln1:ln8, ln9), 2, 0)
## Warning in ln1:ln8: numerical expression has 72 elements: only the first used
## Warning in ln1:ln8: numerical expression has 72 elements: only the first used
nlsite <- rowSums(cbind(ln1, ln2, ln3, ln4, ln5, ln6, ln7, ln8, ln9))
label.var(nlsite, "Number of lymph node sites")
tab1(nlsite, graph=FALSE)
## nlsite : Number of lymph node sites
## Frequency Percent Cum. percent
## 1 17 23.6 23.6
## 2 15 20.8 44.4
## 3 25 34.7 79.2
## 4 8 11.1 90.3
## 5 6 8.3 98.6
## 6 1 1.4 100.0
## Total 72 100.0 100.0
data.frame(case, hn, ln1,ln2, ln3, ln4, ln5, ln6, ln7, ln8, ln9)[nlsite>5,]
## case hn ln1 ln2 ln3 ln4 ln5 ln6 ln7 ln8 ln9
## 30 53 1587899 1 1 1 1 1 0 1 0 0
nlsitegp <- cut(nlsite, br=c(0,2,9), labels=c("<3","3+"))
label.var(nlsitegp, "Number of lymph node sites")
tab1(nlsitegp, graph=FALSE)
## nlsitegp : Number of lymph node sites
## Frequency Percent Cum. percent
## <3 32 44.4 44.4
## 3+ 40 55.6 100.0
## Total 72 100.0 100.0
abdom <- factor(ifelse(ln2==1, "Yes","No"));
label.var(abdom, "Abdominal mass")
mediast <- factor(ifelse(ln1==1, "Yes","No"));
label.var(mediast, "Mediastinal mass")
# Splenomegaly
splenomega <- factor(splenomega, labels=noyes, levels=c(2,1))
label.var(splenomega, "Splenomegaly")
tab1(splenomega, graph=FALSE)
## splenomega : Splenomegaly
## Frequency Percent Cum. percent
## No 57 79.2 79.2
## Yes 15 20.8 100.0
## Total 72 100.0 100.0
summ(if3, graph=FALSE) # Distance below LCM
## obs. mean median s.d. min. max.
## 15 4 3 2.699 1 11
# Hepatomegaly
hepatomega <- factor(hepatomega, labels=noyes, levels=c(2,1))
label.var(hepatomega, "Hepatomegaly")
tab1(hepatomega, graph=FALSE)
## hepatomega : Hepatomegaly
## Frequency Percent Cum. percent
## No 43 59.7 59.7
## Yes 29 40.3 100.0
## Total 72 100.0 100.0
summ(if4, graph=FALSE) # Distance below RCM
## obs. mean median s.d. min. max.
## 29 3.517 3 1.479 1 6
# Staging Work-up
bonc <- factor(bonc, labels=noyes, levels=c(2,1))
label.var(bonc, "Bone scan")
tab1(bonc, graph=FALSE)
## bonc : Bone scan
## Frequency Percent Cum. percent
## No 50 69.4 69.4
## Yes 22 30.6 100.0
## Total 72 100.0 100.0
if5 <- factor(if5, labels=c("Positive","Negative"))
label.var(if5, "Bone scan result")
tab1(if5, graph=FALSE)
## if5 : Bone scan result
## Frequency %(NA+) %(NA-)
## Positive 3 4.2 13.6
## Negative 19 26.4 86.4
## NA's 50 69.4 0.0
## Total 72 100.0 100.0
# All had a bone marrow aspiration
tab1(bma, graph=FALSE)
## bma :
## Frequency Percent Cum. percent
## 1 72 100 100
## Total 72 100 100
if6 <- factor(if6, labels=noyes, levels=c(2,1))
label.var(if6, "Bone marrow involvement")
tab1(if6, graph=FALSE) # Of which 3 had involvement
## if6 : Bone marrow involvement
## Frequency Percent Cum. percent
## No 69 95.8 95.8
## Yes 3 4.2 100.0
## Total 72 100.0 100.0
# Bone or bone marrow involvement
addmargins(table(if5, if6, useNA="ifany"))
## if6
## if5 No Yes Sum
## Positive 2 1 3
## Negative 19 0 19
## <NA> 48 2 50
## Sum 69 3 72
bonemets <- factor(ifelse((bonc=="Yes" & if5=="Positive") | if6=="Yes", "Yes","No"))
label.var(bonemets, "Bone or marrow metastases")
tab1(bonemets, graph=FALSE)
## bonemets : Bone or marrow metastases
## Frequency Percent Cum. percent
## No 67 93.1 93.1
## Yes 5 6.9 100.0
## Total 72 100.0 100.0
ctches <- factor(ctches, labels=noyes, levels=c(2,1))
label.var(ctches, "CT chest")
tab1(ctches, graph=FALSE)
## ctches : CT chest
## Frequency Percent Cum. percent
## No 27 37.5 37.5
## Yes 45 62.5 100.0
## Total 72 100.0 100.0
if7 <- factor(if7, labels=c("Pulmonary mass","Other","Negative"), levels=c(1,3,2))
label.var(if7, "Chest involvement")
tab1(if7, graph=FALSE)
## if7 : Chest involvement
## Frequency %(NA+) %(NA-)
## Pulmonary mass 18 25.0 40.0
## Other 10 13.9 22.2
## Negative 17 23.6 37.8
## NA's 27 37.5 0.0
## Total 72 100.0 100.0
label.var(if8, "Other chest involvement")
tab1(if8, graph=FALSE, sort="decreasing")
## if8 : Other chest involvement
## Frequency %(NA+) %(NA-)
## <NA> 62 86.1 0
## mediastinal mass 5 6.9 50
## axillary LN enlargement 3 4.2 30
## cervical LN enlargement 1 1.4 10
## pericardial effusion, pleural effusion 1 1.4 10
## Total 72 100.0 100
ctab <- factor(ctab, labels=noyes, levels=c(2,1))
label.var(ctab, "CT of abdomen")
tab1(ctab, graph=FALSE)
## ctab : CT of abdomen
## Frequency Percent Cum. percent
## No 12 16.7 16.7
## Yes 60 83.3 100.0
## Total 72 100.0 100.0
if9 <- factor(if9, labels=c("Hepatosplenomegaly","Other","Negative"), levels=c(1,3,2))
label.var(if9, "CT abdomen results")
tab1(if9, graph=FALSE)
## if9 : CT abdomen results
## Frequency %(NA+) %(NA-)
## Hepatosplenomegaly 33 45.8 55.0
## Other 7 9.7 11.7
## Negative 20 27.8 33.3
## NA's 12 16.7 0.0
## Total 72 100.0 100.0
recode.is.na(if10, "_No CT_")
label.var(if10, "Other CT abdomen results")
tab1(if10, graph=FALSE)
## if10 : Other CT abdomen results
## Frequency Percent Cum. percent
## _No CT_ 65 90.3 90.3
## abdominal mass 4 5.6 95.8
## hepatic mass 1 1.4 97.2
## renal nodule 1 1.4 98.6
## splenic nodule 1 1.4 100.0
## Total 72 100.0 100.0
# Other imaging study
other <- factor(other, labels=noyes, levels=c(2,1))
label.var(other, "Other imaging study")
tab1(other, graph=FALSE)
## other : Other imaging study
## Frequency Percent Cum. percent
## No 49 68.1 68.1
## Yes 23 31.9 100.0
## Total 72 100.0 100.0
recode.is.na(if11, "_No Imaging_")
label.var(if11, "Other imaging study results")
tab1(if11, graph=FALSE)
## if11 : Other imaging study results
## Frequency Percent Cum. percent
## _No Imaging_ 49 68.1 68.1
## CT abdomen-abdominal mass 2 2.8 70.8
## CT abdomen-adrenal gland involvement 1 1.4 72.2
## CT abdomen-splenic nodule, adrenal nodule 1 1.4 73.6
## CT neck-cervical LN enlargement 5 6.9 80.6
## CXR-mediastinal mass 1 1.4 81.9
## CXR-normal 2 2.8 84.7
## CXR-perihilar LN enlargement 1 1.4 86.1
## MRI spine-paravertebral mass 1 1.4 87.5
## USG abdomen-abdominal LN enlargement 2 2.8 90.3
## USG abdomen-normal 6 8.3 98.6
## USG abdomen-splenomegaly 1 1.4 100.0
## Total 72 100.0 100.0
# Part III: Investigations ####
# CBC
wbc <- wbc/1000
summ(wbc); abline(v=13.5, lty=2)
## obs. mean median s.d. min. max.
## 72 11.075 10.325 4.4 2.25 23.52
wbc13.5 <- factor(ifelse(wbc>13.5, " > 13.5"," <= 13.5"))
label.var(wbc13.5, "White blood cells")
summ(hb); abline(v=11, lty=2)
## obs. mean median s.d. min. max.
## 72 10.021 10.5 2.184 3.7 14.1
hb11 <- factor(ifelse(hb>11, " > 11.0", " <= 11.0"))
label.var(hb11, "Hemoglobin")
summ(hct)
## obs. mean median s.d. min. max.
## 72 31.418 33 6.013 11 42.2
summ(plt); abline(v=411500, lty=2)
## obs. mean median s.d. min. max.
## 72 435958.333 411500 168153.136 45000 934000
plt.cnt <- factor(ifelse(plt<400000, "< 400,000","400,000+"))
label.var(plt.cnt, "Platelet count")
# Serum chemistry
summ(bun)
## obs. mean median s.d. min. max.
## 72 9.608 9.1 3.12 4.2 18.5
summ(creatinine)
## obs. mean median s.d. min. max.
## 72 0.528 0.515 0.16 0.24 0.99
summ(potassium)
## obs. mean median s.d. min. max.
## 72 4.002 4 0.482 2.41 5.4
summ(bicarbonat)
## obs. mean median s.d. min. max.
## 72 23.008 23 2.409 17 29
summ(calcium)
## obs. mean median s.d. min. max.
## 72 9.276 9.45 0.736 6.8 10.5
summ(phosphorus)
## obs. mean median s.d. min. max.
## 72 4.653 4.7 0.793 2.8 6.6
summ(uric)
## obs. mean median s.d. min. max.
## 71 4.335 4.2 1.711 1.8 13.9
summ(sgot)
## obs. mean median s.d. min. max.
## 72 24.708 24 10.614 3 66
summ(sgpt)
## obs. mean median s.d. min. max.
## 72 18.306 14 12.939 3 62
summ(alp)
## obs. mean median s.d. min. max.
## 72 232.861 206 137.758 35 911
summ(alb)
## obs. mean median s.d. min. max.
## 72 3.911 4.05 0.715 1.9 5.3
alb3.5 <- factor(ifelse(alb<3.5," < 3.5"," 3.5+"))
label.var(alb3.5, " Albumin level")
summ(esr)
## obs. mean median s.d. min. max.
## 6 55.167 50.5 29.24 18 98
ldh <- factor(ldh, labels=noyes)
label.var(ldh, "LDH")
tab1(ldh, graph=FALSE)
## ldh : LDH
## Frequency Percent Cum. percent
## No 62 86.1 86.1
## Yes 10 13.9 100.0
## Total 72 100.0 100.0
summ(if12, main="Distribution of Lactate Dehydrogenase")
## obs. mean median s.d. min. max.
## 62 579.194 501 268.247 306 1657
# Part IV: Diagnosis ####
dx <- factor(dx, labels=c("Classical Hodgkin's'","NLP Hodgkin's"))
label.var(dx, "Histological subtype")
tab1(dx, graph=FALSE)
## dx : Histological subtype
## Frequency Percent Cum. percent
## Classical Hodgkin's' 68 94.4 94.4
## NLP Hodgkin's 4 5.6 100.0
## Total 72 100.0 100.0
chd <- factor(chd, labels=c("Mixed cellularity","lymphocyte depletion","lymphocytic predominance","Nodular sclerosis"))
label.var(chd, "Classical subtype")
tab1(chd, graph=FALSE)
## chd : Classical subtype
## Frequency %(NA+) %(NA-)
## Mixed cellularity 28 38.9 41.2
## lymphocyte depletion 7 9.7 10.3
## lymphocytic predominance 3 4.2 4.4
## Nodular sclerosis 30 41.7 44.1
## NA's 4 5.6 0.0
## Total 72 100.0 100.0
risk <- factor(risk, labels=c("Low","Intermediate","High"))
label.var(risk, "Risk group")
tab1(risk, graph=FALSE)
## risk : Risk group
## Frequency Percent Cum. percent
## Low 12 16.7 16.7
## Intermediate 36 50.0 66.7
## High 24 33.3 100.0
## Total 72 100.0 100.0
ebv <- factor(ebv, labels=noyes, levels=c(2,1))
label.var(ebv, "EBV infection work-up")
tab1(ebv, graph=FALSE)
## ebv : EBV infection work-up
## Frequency Percent Cum. percent
## No 67 93.1 93.1
## Yes 5 6.9 100.0
## Total 72 100.0 100.0
if13 <- factor(if13, labels=c("positive","negative"))
label.var(if13, "EBV infection result")
tab1(if13, graph=FALSE)
## if13 : EBV infection result
## Frequency %(NA+) %(NA-)
## positive 2 2.8 40
## negative 3 4.2 60
## NA's 67 93.1 0
## Total 72 100.0 100
# Part V: Staging ####
stage <- factor(if14, labels=c("Stage I","Stage Ie","Stage II","Stage IIe","Stage III","Stage IIIs","Stage IIIe","Stage IIIse","Stage IV"), levels=1:9)
label.var(stage, "Stage")
tab1(stage, graph=FALSE)
## stage : Stage
## Frequency Percent Cum. percent
## Stage I 4 5.6 5.6
## Stage Ie 0 0.0 5.6
## Stage II 18 25.0 30.6
## Stage IIe 2 2.8 33.3
## Stage III 12 16.7 50.0
## Stage IIIs 15 20.8 70.8
## Stage IIIe 2 2.8 73.6
## Stage IIIse 6 8.3 81.9
## Stage IV 13 18.1 100.0
## Total 72 100.0 100.0
stagegp <- factor(if14, labels=c("Stage I","Stage I","Stage II","Stage II","Stage III","Stage III","Stage III","Stage III","Stage IV"), levels=1:9)
label.var(stagegp, "Stage")
tab1(stagegp, graph=FALSE)
## stagegp : Stage
## Frequency Percent Cum. percent
## Stage I 4 5.6 5.6
## Stage II 20 27.8 33.3
## Stage III 35 48.6 81.9
## Stage IV 13 18.1 100.0
## Total 72 100.0 100.0
bulk <- factor(bulk, labels=noyes, levels=c(2,1))
label.var(bulk, "Bulky disease")
tab1(bulk, graph=FALSE)
## bulk : Bulky disease
## Frequency Percent Cum. percent
## No 51 70.8 70.8
## Yes 21 29.2 100.0
## Total 72 100.0 100.0
# Part VI: Treatment ####
tab1(n1st1, graph=FALSE) # All started first-line chemo
## n1st1 :
## Frequency Percent Cum. percent
## 1 72 100 100
## Total 72 100 100
prot1 <- factor(if15, labels=c("Low risk","Inter/high risk","Other"))
label.var(prot1, "First protocol")
tab1(prot1, graph=FALSE)
## prot1 : First protocol
## Frequency Percent Cum. percent
## Low risk 1 1.4 1.4
## Inter/high risk 8 11.1 12.5
## Other 63 87.5 100.0
## Total 72 100.0 100.0
recode.is.na(if16, "Other")
if16 <- factor(if16, levels=c("COMP","COPP","COPP-ABV","DAL-HD-90","HD-01-05","HD-01-08", "HD-02-05","HD-02-08","HD-03-05","Modified COPP","Modified MOPP","PPCV", "Other"))
tab1(if16, graph=FALSE)
## if16 :
## Frequency Percent Cum. percent
## COMP 1 1.4 1.4
## COPP 3 4.2 5.6
## COPP-ABV 13 18.1 23.6
## DAL-HD-90 2 2.8 26.4
## HD-01-05 1 1.4 27.8
## HD-01-08 3 4.2 31.9
## HD-02-05 1 1.4 33.3
## HD-02-08 18 25.0 58.3
## HD-03-05 5 6.9 65.3
## Modified COPP 1 1.4 66.7
## Modified MOPP 14 19.4 86.1
## PPCV 1 1.4 87.5
## Other 9 12.5 100.0
## Total 72 100.0 100.0
chemo <- if16
label.var(chemo, "Chemotherapy protocol")
recode(chemo, c("HD-01-05","HD-02-05","HD-03-05"), "HD-0X-05")
recode(chemo, c("HD-01-08","HD-02-08","HD-03-08"), "HD-0X-08")
## Warning in recode.default(chemo, c("HD-01-08", "HD-02-08", "HD-03-08"), : The old.value is/are not element of levels of 'chemo'
recode(chemo, c("COPP-ABV","COMP", "Modified COPP"), "COPP")
recode(chemo, c("DAL-HD-90","PPCV"), "Other")
tab1(chemo, graph=FALSE)
## chemo : Chemotherapy protocol
## Frequency Percent Cum. percent
## COPP 18 25.0 25.0
## Other 12 16.7 41.7
## HD-0X-05 7 9.7 51.4
## HD-0X-08 21 29.2 80.6
## Modified MOPP 14 19.4 100.0
## Total 72 100.0 100.0
start <- as.Date(start)
label.var(start, "Date started chemotherapy")
summ(start)
## obs. mean median s.d. min. max.
## 72 2547-09-19 2550-06-01 <NA> 2528-05-13 2560-08-22
era <- factor(ifelse(start<as.Date("2545-01-01"), "2528 - 2544","2545 - 2560"))
label.var(era, "Era of diagnosis")
tab1(era, graph=FALSE)
## era : Era of diagnosis
## Frequency Percent Cum. percent
## 2528 - 2544 23 31.9 31.9
## 2545 - 2560 49 68.1 100.0
## Total 72 100.0 100.0
rt1 <- factor(rt1, labels=noyes, levels=c(2,1))
label.var(rt1, "Radiotherapy")
tab1(rt1, graph=FALSE)
## rt1 : Radiotherapy
## Frequency Percent Cum. percent
## No 27 37.5 37.5
## Yes 45 62.5 100.0
## Total 72 100.0 100.0
complete1 <- factor(complete1, labels=noyes, levels=c(2,1))
label.var(complete1, "Completed protocol")
tab1(complete1, graph=FALSE)
## complete1 : Completed protocol
## Frequency Percent Cum. percent
## No 20 27.8 27.8
## Yes 52 72.2 100.0
## Total 72 100.0 100.0
if17 <- as.Date(if17)
label.var(if17, "Date chemotherapy completed")
summ(if17)
## obs. mean median s.d. min. max.
## 52 2549-06-02 2550-11-21 <NA> 2529-11-17 2560-12-19
# Duration of protocol therapy
prot.time <- as.numeric(difftime(if17, start))/12
label.var(prot.time, "Duration of treatment (months)")
summ(prot.time)
## obs. mean median s.d. min. max.
## 52 21.372 14.583 14.717 7.417 63.917
sortBy(start)
plot(start, 1:72, pch=18)
points(if17, 1:72, pch=18, col="red")
segments(start, 1:72, if17, 1:72)
# Response
respon1 <- factor(respon1, labels=c("Complete","Partial","No change","PD","Death"))
label.var(respon1, "Response to first protocol")
tab1(respon1, graph=FALSE)
## respon1 : Response to first protocol
## Frequency Percent Cum. percent
## Complete 48 66.7 66.7
## Partial 15 20.8 87.5
## No change 2 2.8 90.3
## PD 5 6.9 97.2
## Death 2 2.8 100.0
## Total 72 100.0 100.0
if4pd1dt <- as.Date(if4pd1dt)
label.var(if4pd1dt, "Date of PD")
summ(if4pd1dt, graph=FALSE)
## obs. mean median s.d. min. max.
## 5 2546-10-09 2547-11-26 <NA> 2537-02-22 2556-12-08
relap <- factor(relap, labels=noyes, levels=c(2,1))
label.var(relap, "Relapsed")
tab1(relap, graph=FALSE)
## relap : Relapsed
## Frequency Percent Cum. percent
## No 63 87.5 87.5
## Yes 9 12.5 100.0
## Total 72 100.0 100.0
if18 <- as.Date(if18)
label.var(if18, "Date of relapse")
summ(if18, graph=FALSE)
## obs. mean median s.d. min. max.
## 9 2547-08-25 2549-01-09 <NA> 2533-08-15 2553-06-30
st <- factor(st, labels=c("nodal", "extranodal"))
label.var(st, "Site of relapse")
tab1(st, graph=FALSE)
## st : Site of relapse
## Frequency %(NA+) %(NA-)
## nodal 5 6.9 55.6
## extranodal 4 5.6 44.4
## NA's 63 87.5 0.0
## Total 72 100.0 100.0
xten <- factor(xten, labels=c("localised", "extensive"))
label.var(xten, "Extent of relapse")
tab1(xten, graph=FALSE)
## xten : Extent of relapse
## Frequency %(NA+) %(NA-)
## localised 2 2.8 22.2
## extensive 7 9.7 77.8
## NA's 63 87.5 0.0
## Total 72 100.0 100.0
# Second line chemo
n2nd <- factor(n2nd, labels=noyes, levels=c(2,1))
label.var(n2nd, "Started second-line chemo")
tab1(n2nd, graph=FALSE)
## n2nd : Started second-line chemo
## Frequency Percent Cum. percent
## No 59 81.9 81.9
## Yes 13 18.1 100.0
## Total 72 100.0 100.0
start1 <- as.Date(start1)
label.var(start1, "Started second-line chemotherapy")
summ(start1, graph=FALSE)
## obs. mean median s.d. min. max.
## 13 2547-03-04 2548-11-21 <NA> 2533-08-29 2557-01-13
respon2 <- factor(respon2, labels=c("Complete","Partial","No change","PD","Death"))
label.var(respon2, "Response to second protocol")
tab1(respon2, graph=FALSE)
## respon2 : Response to second protocol
## Frequency %(NA+) %(NA-)
## Complete 1 1.4 7.7
## Partial 2 2.8 15.4
## No change 1 1.4 7.7
## PD 5 6.9 38.5
## Death 4 5.6 30.8
## NA's 59 81.9 0.0
## Total 72 100.0 100.0
# Bone marrow transplant
tab1(bmt, graph=FALSE) # 1 = Yes, 2 = No
## bmt :
## Frequency Percent Cum. percent
## 1 2 2.8 2.8
## 2 70 97.2 100.0
## Total 72 100.0 100.0
tab1(bmt1, graph=FALSE) # 1 = Autologous HSCT, 2 = Allogenic HSCT
## bmt1 :
## Frequency %(NA+) %(NA-)
## 1 1 1.4 50
## 2 1 1.4 50
## <NA> 70 97.2 0
## Total 72 100.0 100
# Part VII: Complications ####
tab1(n2nd1, graph=FALSE)
## n2nd1 :
## Frequency Percent Cum. percent
## 1 1 1.4 1.4
## 2 71 98.6 100.0
## Total 72 100.0 100.0
toxic <- factor(toxic, labels=noyes, levels=c(2,1))
label.var(toxic, "Toxicity")
tab1(toxic, graph=FALSE)
## toxic : Toxicity
## Frequency Percent Cum. percent
## No 63 87.5 87.5
## Yes 9 12.5 100.0
## Total 72 100.0 100.0
tableStack(renal:other1) # 1 = Yes, 2 = No
## 1 2 count mean sd description
## renal 0 9 9 2 0
## liver 4 5 9 1.6 0.5
## neuro 1 8 9 1.9 0.3
## blad 2 7 9 1.8 0.4
## muco 1 8 9 1.9 0.3
## cardio 0 9 9 2 0
## other1 3 6 9 1.7 0.5
## Total score 9 12.8 0.7
## Average score 9 1.8 0.1
tab1(if21, graph=FALSE) # Other toxicity
## if21 :
## Frequency %(NA+) %(NA-)
## Hypertensive encephalopathy 1 1.4 33.3
## restrictive lung disease 1 1.4 33.3
## SIADH 1 1.4 33.3
## <NA> 69 95.8 0.0
## Total 72 100.0 100.0
# Part VIII: Current status ####
status <- factor(status, labels=c("Alive, on chemo","Alive, off chemo","Alive, on other Tx","Lost to FU","Dead"), levels=1:5)
tab1(status, graph=FALSE)
## status :
## Frequency Percent Cum. percent
## Alive, on chemo 0 0.0 0.0
## Alive, off chemo 51 70.8 70.8
## Alive, on other Tx 0 0.0 70.8
## Lost to FU 6 8.3 79.2
## Dead 15 20.8 100.0
## Total 72 100.0 100.0
status <- factor(status) # Omit categories that have no patients
label.var(status, "Final outcome status")
cause <- factor(cause, labels=c("Cancer related","Toxicity","Infection"), levels=1:3)
label.var(cause, "Cause of death")
tab1(cause, graph=FALSE)
## cause : Cause of death
## Frequency %(NA+) %(NA-)
## Cancer related 9 12.5 60
## Toxicity 0 0.0 0
## Infection 6 8.3 40
## NA's 57 79.2 0
## Total 72 100.0 100
if22 <- as.Date(if22)
label.var(if22, "Date of death")
summ(if22)
## obs. mean median s.d. min. max.
## 15 2549-03-14 2551-05-11 <NA> 2531-10-22 2558-05-05
last <- as.Date(last)
label.var(last, "Last known date alive")
summ(last)
## obs. mean median s.d. min. max.
## 72 2553-12-20 2556-12-12 <NA> 2531-10-22 2562-01-20
summ(last, by=status)
## For status = Alive, off chemo
## obs. mean median s.d. min. max.
## 51 2556-12-15 2560-01-23 <NA> 2536-05-10 2562-01-20
##
## For status = Lost to FU
## obs. mean median s.d. min. max.
## 6 2540-07-04 2538-09-22 <NA> 2532-03-19 2558-03-13
##
## For status = Dead
## obs. mean median s.d. min. max.
## 15 2549-03-14 2551-05-11 <NA> 2531-10-22 2558-05-05
data.frame(case, hn, start, if17, respon1, last)[status=="Alive, off chemo" & last<as.Date("2540-01-01"),]
## case hn start if17 respon1 last
## 12 13 434706 2534-05-10 2535-11-02 Complete 2536-05-10
## 13 12 462132 2534-11-22 <NA> Complete 2539-03-04
addmargins(table(respon1, status))
## status
## respon1 Alive, off chemo Lost to FU Dead Sum
## Complete 40 3 5 48
## Partial 9 1 5 15
## No change 1 0 1 2
## PD 1 2 2 5
## Death 0 0 2 2
## Sum 51 6 15 72
# Of the 41 patients who got complete response from first-line treatment, 5 died, ...
tab1(cause[respon1=="Complete"], graph=FALSE)
## cause[respon1 == "Complete"] :
## Frequency %(NA+) %(NA-)
## Cancer related 3 6.2 60
## Toxicity 0 0.0 0
## Infection 2 4.2 40
## NA's 43 89.6 0
## Total 48 100.0 100
# (2 disease-related, 3 due to infection), and 36 are alive and well.
table(relap, status, respon1)
## , , respon1 = Complete
##
## status
## relap Alive, off chemo Lost to FU Dead
## No 39 3 1
## Yes 1 0 4
##
## , , respon1 = Partial
##
## status
## relap Alive, off chemo Lost to FU Dead
## No 8 1 2
## Yes 1 0 3
##
## , , respon1 = No change
##
## status
## relap Alive, off chemo Lost to FU Dead
## No 1 0 1
## Yes 0 0 0
##
## , , respon1 = PD
##
## status
## relap Alive, off chemo Lost to FU Dead
## No 1 2 2
## Yes 0 0 0
##
## , , respon1 = Death
##
## status
## relap Alive, off chemo Lost to FU Dead
## No 0 0 2
## Yes 0 0 0
data.frame(case, start, respon1, relap, n2nd, toxic, n2nd1, if20, status, cause, if22)[respon1=="Complete" & relap=="No" & status=="Dead",]
## case start respon1 relap n2nd toxic n2nd1 if20 status cause if22
## 52 55 2553-02-18 Complete No No No 2 <NA> Dead Infection 2553-09-26
# case 33 had a complete response but had a cancer-related death and has no information about relapse (date, site...) or other malignnacy
data.frame(case, relap, n2nd, respon2, toxic, status)[respon1=="Complete" & relap=="Yes",]
## case relap n2nd respon2 toxic status
## 19 20 Yes Yes Death No Dead
## 27 29 Yes Yes Complete No Alive, off chemo
## 28 37 Yes Yes PD No Dead
## 32 33 Yes No <NA> No Dead
## 40 42 Yes Yes No change Yes Dead
# Of the 4 who relapsed after achieving a complete response, 3 died and 1 is alive and well.
addmargins(table(respon1, status))
## status
## respon1 Alive, off chemo Lost to FU Dead Sum
## Complete 40 3 5 48
## Partial 9 1 5 15
## No change 1 0 1 2
## PD 1 2 2 5
## Death 0 0 2 2
## Sum 51 6 15 72
# Of the 16 patients who got partial response from first-line treatment, 5 died, ...
tab1(cause[respon1=="Partial"], graph=FALSE)
## cause[respon1 == "Partial"] :
## Frequency %(NA+) %(NA-)
## Cancer related 4 26.7 80
## Toxicity 0 0.0 0
## Infection 1 6.7 20
## NA's 10 66.7 0
## Total 15 100.0 100
# (4 disease-related, 1 due to infection), and 11 are alive and well.
# Outcome variable definitions
# Outcomes are disease free survival (DFS) and overall survival (OS).
# DFS is defined as the time from starting treatment [start] to date of progressive disease [if4pddate], relapse [if18], death [if22], or most recent follow-up [last].
pddate <- as.Date(ifelse(respon1=="PD", if4pd1dt,
ifelse(relap=="Yes", if18,
ifelse(!is.na(respon2) & respon2=="PD", if4pd2dt,
ifelse(status=="Dead", if22, last)))), origin="1970-01-01")
label.var(pddate, "Date of PD or last non-PD date")
summ(pddate)
## obs. mean median s.d. min. max.
## 72 2553-08-08 2555-12-03 <NA> 2531-10-22 2562-01-20
dfs <- as.numeric(difftime(pddate, start)/30)
label.var(dfs, "Disease-free survival time")
summ(dfs)
## obs. mean median s.d. min. max.
## 72 71.66 41.55 66.552 1.467 286.4
dfs.status <- as.integer(ifelse(respon1=="PD" | relap=="Yes" | (!is.na(respon2) & respon2=="PD") | (status=="Dead" & cause=="Cancer related"), 1, 0))
label.var(dfs.status, "Disease-free survival status")
tab1(dfs.status, graph=FALSE)
## dfs.status : Disease-free survival status
## Frequency Percent Cum. percent
## 0 55 76.4 76.4
## 1 17 23.6 100.0
## Total 72 100.0 100.0
addmargins(table(status, dfs.status))
## dfs.status
## status 0 1 Sum
## Alive, off chemo 48 3 51
## Lost to FU 4 2 6
## Dead 3 12 15
## Sum 55 17 72
data.frame(case, hn, start, respon1, last, status, cause)[which(status=="Dead" & dfs.status==0),]
## case hn start respon1 last status cause
## 4 5 287053 2531-08-18 Death 2531-10-22 Dead Infection
## 52 55 1617834 2553-02-18 Complete 2553-09-26 Dead Infection
## 63 66 1929445 2557-12-07 Death 2558-05-05 Dead Infection
#install.packages("survminer")
library(survminer)
## Warning: package 'survminer' was built under R version 3.5.2
## Loading required package: ggplot2
##
## Attaching package: 'ggplot2'
## The following object is masked _by_ '.GlobalEnv':
##
## .data
## The following object is masked from 'package:epicalc':
##
## alpha
## Loading required package: ggpubr
## Warning: package 'ggpubr' was built under R version 3.5.2
## Loading required package: magrittr
# Disease-free survival ####
fit <- survfit(Surv(dfs, dfs.status) ~ 1, data=.data)
summary(fit, times=60) # 5-year disease-free survival rate (70.7%, [59.7% - 83.8%])
## Call: survfit(formula = Surv(dfs, dfs.status) ~ 1, data = .data)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 60 31 17 0.707 0.0611 0.597 0.838
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.13,
tables.theme = theme_cleantable(),
# Color palette
palette = c("#E7B800", "#2E9FDF"),
ggtheme = theme_bw(),
legend="none",
title="Disease-free survival",
xlab="Time (months)"
)
# Disease-free survival rates at 1-year time points.
summary(fit, times=seq(0,60,12))
## Call: survfit(formula = Surv(dfs, dfs.status) ~ 1, data = .data)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 72 0 1.000 0.0000 1.000 1.000
## 12 64 4 0.942 0.0281 0.889 0.999
## 24 50 5 0.864 0.0424 0.784 0.951
## 36 39 5 0.770 0.0548 0.669 0.885
## 48 35 2 0.729 0.0591 0.622 0.854
## 60 31 1 0.707 0.0611 0.597 0.838
# Overall survival ####
os <- as.numeric(difftime(last, start)/30)
label.var(os, "Overall survival time")
summ(os, graph=FALSE)
## obs. mean median s.d. min. max.
## 72 76.128 53.75 65.784 1.467 286.4
os.status <- status=="Dead"
label.var(os.status, "Dead")
tab1(os.status, graph=FALSE)
## os.status : Dead
## Frequency Percent Cum. percent
## FALSE 57 79.2 79.2
## TRUE 15 20.8 100.0
## Total 72 100.0 100.0
summ(os, by=os.status)
## For os.status = FALSE
## obs. mean median s.d. min. max.
## 57 88.065 72 68.421 1.467 286.4
##
## For os.status = TRUE
## obs. mean median s.d. min. max.
## 15 30.767 31.033 20.892 2.167 73.2
addmargins(table(os.status, dfs.status))
## dfs.status
## os.status 0 1 Sum
## FALSE 52 5 57
## TRUE 3 12 15
## Sum 55 17 72
data.frame(case, hn, start, respon1, respon2, last, status, cause)[which(!os.status & dfs.status==1),]
## case hn start respon1 respon2 last status cause
## 11 11 433150 2534-05-01 PD PD 2538-12-12 Lost to FU <NA>
## 18 19 763156 2540-01-30 PD Partial 2540-07-14 Lost to FU <NA>
## 27 29 1197689 2547-04-17 Complete Complete 2561-03-14 Alive, off chemo <NA>
## 41 43 1444174 2550-11-01 Partial Partial 2553-05-12 Alive, off chemo <NA>
## 59 63 1803855 2555-12-17 PD PD 2557-11-22 Alive, off chemo <NA>
fit2 <- survfit(Surv(os, os.status) ~ 1, data = .data)
summary(fit2, times=60)
## Call: survfit(formula = Surv(os, os.status) ~ 1, data = .data)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 60 33 14 0.755 0.0584 0.648 0.878
ggsurvplot(fit2, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.13,
tables.theme = theme_cleantable(),
# Color palette
palette = c("#E7B800", "#2E9FDF"),
ggtheme = theme_bw(),
legend="none",
title="Overall survival"
)
# Survival rates at 1-year time points.
summary(fit2, times=seq(0,60,12))
## Call: survfit(formula = Surv(os, os.status) ~ 1, data = .data)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 72 0 1.000 0.0000 1.000 1.000
## 12 66 4 0.943 0.0276 0.891 0.999
## 24 54 2 0.913 0.0340 0.849 0.982
## 36 45 2 0.877 0.0411 0.800 0.961
## 48 38 5 0.775 0.0561 0.673 0.894
## 60 33 1 0.755 0.0584 0.648 0.878
# 5-year survival rate = 75.5%
source("surv.table.R")
# Table 1 ####
Table1a <- tableStack(vars=c(agem, sex, stagegp, risk, dx, bonemets, bsympt, nlsitegp, bulk, hepatomega, splenomega, abdom, mediast, wbc13.5, hb11, plt.cnt, ldh, alb3.5, rt1, era, chemo), by="none")
Table1a
## Total
## Total 72
##
## Median age
## Old 38 (52.8)
## Young 34 (47.2)
##
## Sex
## Male 57 (79.2)
## Female 15 (20.8)
##
## Stage
## Stage I 4 (5.6)
## Stage II 20 (27.8)
## Stage III 35 (48.6)
## Stage IV 13 (18.1)
##
## Risk group
## Low 12 (16.7)
## Intermediate 36 (50)
## High 24 (33.3)
##
## Histological subtype
## Classical Hodgkin's' 68 (94.4)
## NLP Hodgkin's 4 (5.6)
##
## Bone or marrow metastases
## No 67 (93.1)
## Yes 5 (6.9)
##
## B symptoms
## No 40 (55.6)
## Yes 32 (44.4)
##
## Number of lymph node sites
## <3 32 (44.4)
## 3+ 40 (55.6)
##
## Bulky disease
## No 51 (70.8)
## Yes 21 (29.2)
##
## Hepatomegaly
## No 43 (59.7)
## Yes 29 (40.3)
##
## Splenomegaly
## No 57 (79.2)
## Yes 15 (20.8)
##
## Abdominal mass
## No 40 (55.6)
## Yes 32 (44.4)
##
## Mediastinal mass
## No 53 (73.6)
## Yes 19 (26.4)
##
## White blood cells
## <= 13.5 56 (77.8)
## > 13.5 16 (22.2)
##
## Hemoglobin
## <= 11.0 44 (61.1)
## > 11.0 28 (38.9)
##
## Platelet count
## < 400,000 34 (47.2)
## 400,000+ 38 (52.8)
##
## LDH
## No 62 (86.1)
## Yes 10 (13.9)
##
## Albumin level
## < 3.5 19 (26.4)
## 3.5+ 53 (73.6)
##
## Radiotherapy
## No 27 (37.5)
## Yes 45 (62.5)
##
## Era of diagnosis
## 2528 - 2544 23 (31.9)
## 2545 - 2560 49 (68.1)
##
## Chemotherapy protocol
## COPP 18 (25)
## Other 12 (16.7)
## HD-0X-05 7 (9.7)
## HD-0X-08 21 (29.2)
## Modified MOPP 14 (19.4)
##
Table1b <- surv.table(c(agem, sex, stagegp, risk, dx, bonemets, bsympt, nlsitegp, bulk, hepatomega, splenomega, abdom, mediast, wbc13.5, hb11, plt.cnt, ldh, alb3.5, rt1, era, chemo), 60)
print.noquote(Table1b)
## DFS Rate lower upper P-value OS Rate lower upper P-value
## Median age 0.568 0.345
## Old 71.2 55.8 90.7 80.7 66.6 97.7
## Young 69.8 55.0 88.7 70.0 55.1 88.9
##
## Sex 0.063 0.284
## Male 64.5 51.5 80.7 72.0 59.4 87.3
## Female 91.7 77.3 100.0 85.6 68.8 100.0
##
## Stage 0.126 0.476
## Stage I 100.0 100.0 100.0 100.0 100.0 100.0
## Stage II 86.9 71.4 100.0 83.1 67.2 100.0
## Stage III 58.6 43.4 79.2 70.9 55.6 90.4
## Stage IV 68.8 38.1 100.0 63.5 34.4 100.0
##
## Risk group 0 0.108
## Low 100.0 100.0 100.0 90.9 75.4 100.0
## Intermediate 79.8 65.2 97.8 82.1 68.6 98.2
## High 43.3 26.2 71.6 55.4 36.2 84.8
##
## Histological subtype 0.31 0.364
## Classical Hodgkin's' 69.2 57.8 82.9 74.3 63.3 87.2
## NLP Hodgkin's 100.0 100.0 100.0 100.0 100.0 100.0
##
## Bone or marrow metastases 0.324 0.39
## No 69.5 58.2 82.9 74.5 63.6 87.2
## Yes 100.0 100.0 100.0 100.0 100.0 100.0
##
## B symptoms 0.002 0.135
## No 85.9 73.7 100.0 83.3 70.7 98.2
## Yes 52.5 36.5 75.5 65.3 49.0 87.0
##
## Number of lymph node sites 0.901 0.926
## <3 70.5 54.9 90.5 76.1 61.0 95.0
## 3+ 71.2 56.7 89.4 75.4 61.6 92.4
##
## Bulky disease 0.821 0.915
## No 73.9 61.6 88.7 76.6 64.6 90.7
## Yes 63.9 44.2 92.5 73.3 54.0 99.5
##
## Hepatomegaly 0.681 0.605
## No 72.3 58.9 88.7 74.8 61.7 90.7
## Yes 68.4 50.7 92.2 78.1 62.1 98.1
##
## Splenomegaly 0 0.054
## No 80.6 69.9 93.0 78.2 67.0 91.4
## Yes 31.3 12.6 77.2 62.7 38.6 100.0
##
## Abdominal mass 0.504 0.862
## No 74.4 60.2 91.9 78.0 64.5 94.2
## Yes 66.4 50.4 87.4 72.5 56.8 92.5
##
## Mediastinal mass 0.992 0.921
## No 71.5 58.8 87.0 76.5 64.5 90.8
## Yes 69.6 50.3 96.2 73.2 53.7 99.9
##
## White blood cells 0.828 0.599
## <= 13.5 72.0 59.9 86.5 78.3 67.0 91.5
## > 13.5 66.5 44.2 100.0 65.3 42.6 100.0
##
## Hemoglobin 0.125 0.875
## <= 11.0 64.1 49.9 82.4 77.1 64.1 92.8
## > 11.0 81.6 66.6 99.9 73.3 56.9 94.5
##
## Platelet count 0.032 0.204
## < 400,000 86.5 74.8 99.9 83.4 70.8 98.1
## 400,000+ 57.3 41.8 78.5 67.8 52.2 88.0
##
## LDH 0.328 0.913
## No 67.9 55.9 82.6 74.5 62.8 88.4
## Yes 88.9 70.6 100.0 80.0 58.7 100.0
##
## Albumin level 0.009 0.013
## < 3.5 50.0 30.1 83.2 53.8 32.7 88.5
## 3.5+ 78.2 66.4 92.1 82.9 72.0 95.5
##
## Radiotherapy 0.3 0.117
## No 60.2 41.8 86.8 68.5 51.4 91.4
## Yes 76.6 64.1 91.5 79.1 66.3 94.5
##
## Era of diagnosis 0.809 0.647
## 2528 - 2544 70.8 51.8 96.9 84.6 69.6 100.0
## 2545 - 2560 70.4 57.4 86.3 71.6 58.5 87.7
##
## Chemotherapy protocol 0.123 0.505
## COPP 59.9 38.6 93.1 57.9 35.1 95.7
## Other 59.3 25.8 100.0 66.7 30.0 100.0
## HD-0X-05 42.9 18.2 100.0 50.0 22.5 100.0
## HD-0X-08 88.5 74.8 100.0 84.4 69.5 100.0
## Modified MOPP 74.6 53.3 100.0 85.7 69.2 100.0
##
write.csv(Table1a, file="Table 1a.csv")
write.csv(Table1b, file="Table 1b.csv")
# Disease-free survival by...
# Sex
survdiff(Surv(dfs, dfs.status) ~ sex)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ sex)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=Male 57 16 12.66 0.879 3.47
## sex=Female 15 1 4.34 2.566 3.47
##
## Chisq= 3.5 on 1 degrees of freedom, p= 0.06
fit.sex <- survfit(Surv(dfs, dfs.status) ~ sex, data = .data)
ggsurvplot(fit.sex, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.16,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by sex"
)
# Stage
survdiff(Surv(dfs, dfs.status) ~ stagegp)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ stagegp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stagegp=Stage I 4 0 1.12 1.1186 1.199
## stagegp=Stage II 20 2 5.24 2.0007 2.912
## stagegp=Stage III 35 13 8.38 2.5458 5.026
## stagegp=Stage IV 13 2 2.26 0.0307 0.036
##
## Chisq= 5.7 on 3 degrees of freedom, p= 0.1
fit <- survfit(Surv(dfs, dfs.status) ~ stagegp, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by stage",
legend.labs=levels(stagegp)
)
# Risk group
survdiff(Surv(dfs, dfs.status) ~ risk)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ risk)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## risk=Low 12 0 3.34 3.34 4.21
## risk=Intermediate 36 5 8.86 1.68 3.52
## risk=High 24 12 4.80 10.82 15.16
##
## Chisq= 16 on 2 degrees of freedom, p= 0.0003
fit <- survfit(Surv(dfs, dfs.status) ~ risk, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.theme = theme_cleantable(),
# Color palette
palette = c("green3","steelblue","red"),
ggtheme = theme_bw(),
title="Disease-free survival by risk level",
legend.labs=levels(risk)
)
survdiff(Surv(dfs, dfs.status) ~ dx)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ dx)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## dx=Classical Hodgkin's' 68 17 16.029 0.0588 1.03
## dx=NLP Hodgkin's 4 0 0.971 0.9708 1.03
##
## Chisq= 1 on 1 degrees of freedom, p= 0.3
survdiff(Surv(dfs, dfs.status) ~ nlsitegp)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ nlsitegp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## nlsitegp=<3 32 8 7.74 0.00847 0.0156
## nlsitegp=3+ 40 9 9.26 0.00708 0.0156
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9
survdiff(Surv(dfs, dfs.status) ~ bonemets)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ bonemets)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## bonemets=No 67 17 16.087 0.0518 0.974
## bonemets=Yes 5 0 0.913 0.9131 0.974
##
## Chisq= 1 on 1 degrees of freedom, p= 0.3
survdiff(Surv(dfs, dfs.status) ~ bsympt)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ bsympt)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## bsympt=No 40 4 10.14 3.72 9.25
## bsympt=Yes 32 13 6.86 5.49 9.25
##
## Chisq= 9.2 on 1 degrees of freedom, p= 0.002
fit <- survfit(Surv(dfs, dfs.status) ~ bsympt, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.16,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by B symptoms"
)
survdiff(Surv(dfs, dfs.status) ~ bulk)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ bulk)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## bulk=No 51 11 11.44 0.0167 0.051
## bulk=Yes 21 6 5.56 0.0342 0.051
##
## Chisq= 0.1 on 1 degrees of freedom, p= 0.8
survdiff(Surv(dfs, dfs.status) ~ hepatomega)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ hepatomega)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## hepatomega=No 43 10 10.82 0.0615 0.169
## hepatomega=Yes 29 7 6.18 0.1075 0.169
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.7
survdiff(Surv(dfs, dfs.status) ~ splenomega)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ splenomega)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## splenomega=No 57 9 14.38 2.01 13.1
## splenomega=Yes 15 8 2.62 11.06 13.1
##
## Chisq= 13.1 on 1 degrees of freedom, p= 0.0003
fit <- survfit(Surv(dfs, dfs.status) ~ splenomega, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.16,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by splenomegaly"
)
survdiff(Surv(dfs, dfs.status) ~ ln2) # abdominal mass
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ ln2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ln2=0 40 8 9.37 0.200 0.446
## ln2=1 32 9 7.63 0.246 0.446
##
## Chisq= 0.4 on 1 degrees of freedom, p= 0.5
survdiff(Surv(dfs, dfs.status) ~ ln1) # mediastinal mass
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ ln1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ln1=0 53 12 11.98 0.0000319 0.000108
## ln1=1 19 5 5.02 0.0000762 0.000108
##
## Chisq= 0 on 1 degrees of freedom, p= 1
survdiff(Surv(dfs, dfs.status) ~ wbc13.5)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ wbc13.5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## wbc13.5= <= 13.5 56 13 13.37 0.0101 0.0474
## wbc13.5= > 13.5 16 4 3.63 0.0372 0.0474
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8
survdiff(Surv(dfs, dfs.status) ~ hb11)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ hb11)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## hb11= <= 11.0 44 13 9.88 0.982 2.35
## hb11= > 11.0 28 4 7.12 1.364 2.35
##
## Chisq= 2.4 on 1 degrees of freedom, p= 0.1
survdiff(Surv(dfs, dfs.status) ~ plt.cnt)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ plt.cnt)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## plt.cnt=< 400,000 34 4 8.42 2.32 4.6
## plt.cnt=400,000+ 38 13 8.58 2.27 4.6
##
## Chisq= 4.6 on 1 degrees of freedom, p= 0.03
fit <- survfit(Surv(dfs, dfs.status) ~ plt.cnt, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.16,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by platelet count"
)
survdiff(Surv(dfs, dfs.status) ~ ldh)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ ldh)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ldh=No 62 16 14.6 0.135 0.955
## ldh=Yes 10 1 2.4 0.819 0.955
##
## Chisq= 1 on 1 degrees of freedom, p= 0.3
survdiff(Surv(dfs, dfs.status) ~ alb3.5)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ alb3.5)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## alb3.5= < 3.5 19 8 3.6 5.39 6.9
## alb3.5= 3.5+ 53 9 13.4 1.45 6.9
##
## Chisq= 6.9 on 1 degrees of freedom, p= 0.009
fit <- survfit(Surv(dfs, dfs.status) ~ alb3.5, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.16,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by albumin level"
)
survdiff(Surv(dfs, dfs.status) ~ rt1)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ rt1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rt1=No 27 8 5.96 0.695 1.07
## rt1=Yes 45 9 11.04 0.376 1.07
##
## Chisq= 1.1 on 1 degrees of freedom, p= 0.3
# Cox regression for chemistry values
cox.display(coxph(Surv(dfs, dfs.status) ~ bun, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## bun (cont. var.) 1.02 (0.87,1.2) 0.809 0.81
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ creatinine, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## creatinine (cont. var.) 0.81 (0.04,16.52) 0.888 0.888
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ potassium, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## potassium (cont. var.) 1.78 (0.67,4.7) 0.244 0.245
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ bicarbonat, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## bicarbonat (cont. var.) 1.03 (0.84,1.27) 0.767 0.767
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ calcium, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## calcium (cont. var.) 0.92 (0.48,1.76) 0.806 0.807
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ phosphorus, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## phosphorus (cont. var.) 1.52 (0.8,2.89) 0.2 0.2
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ uric, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## uric (cont. var.) 0.83 (0.59,1.17) 0.285 0.242
##
## No. of observations = 71
cox.display(coxph(Surv(dfs, dfs.status) ~ sgot, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## sgot (cont. var.) 1.03 (0.98,1.07) 0.268 0.28
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ sgpt, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## sgpt (cont. var.) 1.03 (0.99,1.06) 0.104 0.127
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ alp, data=.data))
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## alp (cont. var.) 1.001 (0.9982,1.0039) 0.487 0.517
##
## No. of observations = 72
cox.display(coxph(Surv(dfs, dfs.status) ~ stagegp, data=.data))
## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, : Loglik converged before variable 1,2,3 ; beta may be infinite.
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## Stage: ref.=Stage I 0.069
## Stage II 29739026.46 (0,Inf) 0.998
## Stage III 121308413.72 (0,Inf) 0.998
## Stage IV 69599387.28 (0,Inf) 0.998
##
## No. of observations = 72
# Check for colinearity
table(stagegp, risk)
## risk
## stagegp Low Intermediate High
## Stage I 3 1 0
## Stage II 9 11 0
## Stage III 0 17 18
## Stage IV 0 7 6
chisq.test(stagegp, risk)
## Warning in chisq.test(stagegp, risk): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: stagegp and risk
## X-squared = 38.532, df = 6, p-value = 0.000000884
# Multivariate analysis ####
risk2 <- risk; label.var(risk2, "Risk group")
recode(risk2, c("Low","Intermediate"), "Low/inter")
cox0 <- coxph(Surv(dfs, dfs.status) ~ agem+sex+risk2+bulk+bsympt+splenomega+plt.cnt+alb3.5, data=.data)
cox.display(cox0)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Median age: Young vs Old 1.32 (0.51,3.42) 1.67 (0.62,4.5) 0.314 0.313
##
## Sex: Female vs Male 0.18 (0.02,1.37) 0.23 (0.03,1.83) 0.166 0.095
##
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 7.22 (0.95,54.88) 0.056 0.039
##
## Bulky disease: Yes vs No 1.12 (0.41,3.04) 0.45 (0.14,1.41) 0.172 0.158
##
## B symptoms: Yes vs No 4.83 (1.57,14.84) 1.64 (0.22,12.09) 0.63 0.639
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 2.96 (0.91,9.63) 0.072 0.073
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 7.86 (1.78,34.76) 0.007 0.003
##
## Albumin level: 3.5+ vs < 3.5 0.3 (0.11,0.78) 2.46 (0.6,10.09) 0.211 0.209
##
## No. of observations = 72
# Omit age
cox1 <- coxph(Surv(dfs, dfs.status) ~ sex+risk2+bulk+bsympt+splenomega+plt.cnt+alb3.5, data=.data)
cox.display(cox1)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Sex: Female vs Male 0.18 (0.02,1.37) 0.23 (0.03,1.88) 0.171 0.099
##
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 7.53 (1.08,52.41) 0.041 0.025
##
## Bulky disease: Yes vs No 1.12 (0.41,3.04) 0.45 (0.14,1.4) 0.169 0.156
##
## B symptoms: Yes vs No 4.83 (1.57,14.84) 1.38 (0.2,9.48) 0.741 0.746
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 2.96 (0.91,9.58) 0.07 0.071
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 6.97 (1.72,28.29) 0.007 0.003
##
## Albumin level: 3.5+ vs < 3.5 0.3 (0.11,0.78) 2.31 (0.58,9.14) 0.234 0.234
##
## No. of observations = 72
# Omit B symptoms (highly confounded by risk gorup)
table(bsympt, risk2)
## risk2
## bsympt Low/inter High
## No 38 2
## Yes 10 22
cox1a <- coxph(Surv(dfs, dfs.status) ~ sex+risk2+bulk+ splenomega+plt.cnt+alb3.5, data=.data)
cox.display(cox1a)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Sex: Female vs Male 0.18 (0.02,1.37) 0.22 (0.03,1.79) 0.158 0.088
##
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 9.56 (2.42,37.69) 0.001 0.001
##
## Bulky disease: Yes vs No 1.12 (0.41,3.04) 0.44 (0.14,1.37) 0.158 0.145
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 2.9 (0.91,9.31) 0.073 0.074
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 6.9 (1.72,27.71) 0.007 0.003
##
## Albumin level: 3.5+ vs < 3.5 0.3 (0.11,0.78) 2.21 (0.57,8.61) 0.253 0.252
##
## No. of observations = 72
# Omit sex
cox1b <- coxph(Surv(dfs, dfs.status) ~ risk2+bulk+ splenomega+plt.cnt+alb3.5, data=.data)
cox.display(cox1b)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 9.05 (2.28,35.88) 0.002 0.001
##
## Bulky disease: Yes vs No 1.12 (0.41,3.04) 0.38 (0.12,1.18) 0.095 0.085
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 3.79 (1.16,12.33) 0.027 0.029
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 7.01 (1.74,28.32) 0.006 0.003
##
## Albumin level: 3.5+ vs < 3.5 0.3 (0.11,0.78) 2.03 (0.53,7.78) 0.3 0.298
##
## No. of observations = 72
# Omit albumin level
cox1c <- coxph(Surv(dfs, dfs.status) ~ sex+risk2+bulk+ splenomega+plt.cnt, data=.data)
cox.display(cox1c)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Sex: Female vs Male 0.18 (0.02,1.37) 0.24 (0.03,1.89) 0.175 0.102
##
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 6.59 (1.87,23.23) 0.003 0.002
##
## Bulky disease: Yes vs No 1.12 (0.41,3.04) 0.4 (0.12,1.28) 0.121 0.11
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 2.32 (0.77,7.01) 0.134 0.136
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 4.82 (1.4,16.52) 0.012 0.007
##
## No. of observations = 72
# Omit sex and albumin
cox1d <- coxph(Surv(dfs, dfs.status) ~ risk2+bulk+ splenomega+plt.cnt, data=.data)
cox.display(cox1d)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 6.66 (1.84,24.14) 0.004 0.002
##
## Bulky disease: Yes vs No 1.12 (0.41,3.04) 0.36 (0.11,1.17) 0.09 0.08
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 2.98 (1,8.9) 0.051 0.054
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 5.01 (1.46,17.2) 0.01 0.005
##
## No. of observations = 72
# Buly disease ?
table(bulk, risk2)
## risk2
## bulk Low/inter High
## No 39 12
## Yes 9 12
fisher.test(bulk, risk2)
##
## Fisher's Exact Test for Count Data
##
## data: bulk and risk2
## p-value = 0.01197
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.290325 14.648369
## sample estimates:
## odds ratio
## 4.23276
# Omit bulk
cox1e <- coxph(Surv(dfs, dfs.status) ~ risk2+ splenomega+plt.cnt, data=.data)
cox.display(cox1e)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 4.17 (1.37,12.66) 0.012 0.009
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 3.26 (1.16,9.13) 0.025 0.028
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 3.55 (1.13,11.13) 0.03 0.019
##
## No. of observations = 72
# Add stage group (combine stage I + II, and stage III + IV)
stagegp2 <- stagegp; label.var(stagegp2, "Stage group")
recode(stagegp2, "Stage I", "Stage II")
recode(stagegp2, "Stage III", "Stage IV")
cox1f <- coxph(Surv(dfs, dfs.status) ~stagegp2 + risk2+ splenomega+plt.cnt, data=.data)
cox.display(cox1f)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## crude HR(95%CI) adj. HR(95%CI) P(Wald's test) P(LR-test)
## Stage group: Stage IV vs Stage II 4.52 (1.03,19.83) 1.4 (0.22,8.81) 0.719 0.717
##
## Risk group: High vs Low/inter 6.18 (2.17,17.6) 3.65 (1.01,13.14) 0.048 0.028
##
## Splenomegaly: Yes vs No 4.94 (1.9,12.86) 3.12 (1.1,8.85) 0.033 0.035
##
## Platelet count: 400,000+ vs < 400,000 3.19 (1.04,9.8) 3.55 (1.13,11.13) 0.03 0.018
##
## No. of observations = 72
# Final model
model.final <- cox.display(cox1e, crude=FALSE)
write.csv(model.final$table, file="Final.model.csv")
# Table 3 ####
# Model 1 ####
score1 <- rowSums(cbind( (sex=="Male"), (stagegp %in% c("Stage III", "Stage IV")), (bulk=="Yes"), (wbc>13.5), (hb<11.0) ))
label.var(score1, "Model 1 score")
tab1(score1, graph=FALSE)
## score1 : Model 1 score
## Frequency Percent Cum. percent
## 0 2 2.8 2.8
## 1 12 16.7 19.4
## 2 20 27.8 47.2
## 3 26 36.1 83.3
## 4 8 11.1 94.4
## 5 4 5.6 100.0
## Total 72 100.0 100.0
# Combine 0 with 1 and 5 with 4 due to small sample sizes
recode(score1, 0, 1)
recode(score1, 5, 4)
tab1(score1, graph=FALSE)
## score1 : Model 1 score
## Frequency Percent Cum. percent
## 1 14 19.4 19.4
## 2 20 27.8 47.2
## 3 26 36.1 83.3
## 4 12 16.7 100.0
## Total 72 100.0 100.0
survdiff(Surv(dfs, dfs.status) ~ score1)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ score1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## score1=1 14 0 4.12 4.119 5.496
## score1=2 20 3 4.85 0.705 0.988
## score1=3 26 10 5.16 4.543 6.559
## score1=4 12 4 2.87 0.442 0.533
##
## Chisq= 9.9 on 3 degrees of freedom, p= 0.02
fit <- survfit(Surv(dfs, dfs.status) ~ score1, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by Model 1 score"
)
# Combine 3 and 4 ?
#recode(score1, 4, 3)
#tab1(score1, graph=FALSE)
# By response to treatment
respon3 <- respon1; label.var(respon3, "Response to treatment")
recode(respon3, c("No change", "Death"), "PD")
tabpct(score1, respon3, graph=FALSE)
##
## Original table
## Response to treatment
## Model 1 score Complete Partial PD Total
## 1 12 1 1 14
## 2 14 4 2 20
## 3 16 5 5 26
## 4 6 5 1 12
## Total 48 15 9 72
##
## Row percent
## Response to treatment
## Model 1 score Complete Partial PD Total
## 1 12 1 1 14
## (85.7) (7.1) (7.1) (100)
## 2 14 4 2 20
## (70) (20) (10) (100)
## 3 16 5 5 26
## (61.5) (19.2) (19.2) (100)
## 4 6 5 1 12
## (50) (41.7) (8.3) (100)
##
## Column percent
## Response to treatment
## Model 1 score Complete % Partial % PD %
## 1 12 (25.0) 1 (6.7) 1 (11.1)
## 2 14 (29.2) 4 (26.7) 2 (22.2)
## 3 16 (33.3) 5 (33.3) 5 (55.6)
## 4 6 (12.5) 5 (33.3) 1 (11.1)
## Total 48 (100) 15 (100) 9 (100)
score1 <- factor(score1); pack()
surv.table(score1, 60)
## DFS Rate lower upper P-value OS Rate lower upper P-value
## Model 1 score "" "" "" "0.019" "" "" "" "0.386"
## 1 "100.0" "100.0" "100.0" "" "92.9" "80.3" "100.0" ""
## 2 "78.4" "59.3" "100.0" "" "75.0" "56.0" "100.0" ""
## 3 "54.7" "37.0" "81.0" "" "63.6" "44.9" "90.1" ""
## 4 "54.3" "28.8" "100.0" "" "71.4" "44.7" "100.0" ""
## "" "" "" "" "" "" "" ""
# CHIPS ####
score2 <- rowSums(cbind( (stagegp == "Stage IV"), (mediast == "Yes"), (alb<3.5), (fever=="Yes") ))
label.var(score2, "CHIPS score")
tab1(score2, graph=FALSE)
## score2 : CHIPS score
## Frequency Percent Cum. percent
## 0 28 38.9 38.9
## 1 19 26.4 65.3
## 2 18 25.0 90.3
## 3 7 9.7 100.0
## Total 72 100.0 100.0
survdiff(Surv(dfs, dfs.status) ~ score2)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ score2)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## score2=0 28 4 7.49 1.62472 2.92776
## score2=1 19 4 4.17 0.00706 0.00937
## score2=2 18 6 3.85 1.20109 1.55555
## score2=3 7 3 1.49 1.52790 1.68010
##
## Chisq= 4.4 on 3 degrees of freedom, p= 0.2
fit <- survfit(Surv(dfs, dfs.status) ~ score2, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by CHIPS score"
)
# By response to treatment
tab1(score2)
## score2 : CHIPS score
## Frequency Percent Cum. percent
## 0 28 38.9 38.9
## 1 19 26.4 65.3
## 2 18 25.0 90.3
## 3 7 9.7 100.0
## Total 72 100.0 100.0
tabpct(score2, respon3, graph=FALSE)
##
## Original table
## Response to treatment
## CHIPS score Complete Partial PD Total
## 0 22 4 2 28
## 1 12 5 2 19
## 2 9 4 5 18
## 3 5 2 0 7
## Total 48 15 9 72
##
## Row percent
## Response to treatment
## CHIPS score Complete Partial PD Total
## 0 22 4 2 28
## (78.6) (14.3) (7.1) (100)
## 1 12 5 2 19
## (63.2) (26.3) (10.5) (100)
## 2 9 4 5 18
## (50) (22.2) (27.8) (100)
## 3 5 2 0 7
## (71.4) (28.6) (0) (100)
##
## Column percent
## Response to treatment
## CHIPS score Complete % Partial % PD %
## 0 22 (45.8) 4 (26.7) 2 (22.2)
## 1 12 (25.0) 5 (33.3) 2 (22.2)
## 2 9 (18.8) 4 (26.7) 5 (55.6)
## 3 5 (10.4) 2 (13.3) 0 (0.0)
## Total 48 (100) 15 (100) 9 (100)
score2 <- factor(score2); pack()
surv.table(score2, 60)
## DFS Rate lower upper P-value OS Rate lower upper P-value
## CHIPS score "" "" "" "0.222" "" "" "" "0.185"
## 0 "81.3" "66.2" "99.9" "" "81.6" "66.6" "99.9" ""
## 1 "72.0" "51.7" "100.0" "" "86.6" "70.6" "100.0" ""
## 2 "61.0" "40.6" "91.7" "" "65.6" "44.7" "96.4" ""
## 3 "45.7" "17.5" "100.0" "" "45.7" "17.5" "100.0" ""
## "" "" "" "" "" "" "" ""
#
cox.model.1 <- coxph(Surv(dfs, dfs.status) ~ score2, data=.data)
cox.display(cox.model.1)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## CHIPS score: ref.=0 0.241
## 1 1.81 (0.45,7.24) 0.403
## 2 2.94 (0.83,10.44) 0.096
## 3 3.81 (0.85,17.11) 0.081
##
## No. of observations = 72
score2f <- factor(score2); pack()
cox.model.2f <- coxph(Surv(dfs, dfs.status) ~ score2f, data=.data)
cox.display(cox.model.2f)
## Cox's proportional hazard model on time ('dfs') to event ('dfs.status')
##
## HR(95%CI) P(Wald's test) P(LR-test)
## score2f: ref.=0 0.241
## 1 1.81 (0.45,7.24) 0.403
## 2 2.94 (0.83,10.44) 0.096
## 3 3.81 (0.85,17.11) 0.081
##
## No. of observations = 72
# Patients with a score of 3 or more have a significantly higher hazard rate of disease progression compared to patients with a score of 0 or 1.
# Model 3 (this data)
score3 <- rowSums(cbind( (risk2=="High"), (splenomega=="Yes"), (plt.cnt=="400,000+") ))
label.var(score3, "Model 3 score")
tab1(score3, graph=FALSE)
## score3 : Model 3 score
## Frequency Percent Cum. percent
## 0 24 33.3 33.3
## 1 25 34.7 68.1
## 2 17 23.6 91.7
## 3 6 8.3 100.0
## Total 72 100.0 100.0
recode(score3, 3, 2) # Combine 3 with 2
tab1(score3, graph=FALSE)
## score3 : Model 3 score
## Frequency Percent Cum. percent
## 0 24 33.3 33.3
## 1 25 34.7 68.1
## 2 23 31.9 100.0
## Total 72 100.0 100.0
survdiff(Surv(dfs, dfs.status) ~ score3)
## Call:
## survdiff(formula = Surv(dfs, dfs.status) ~ score3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## score3=0 24 1 5.96 4.13 6.36
## score3=1 25 4 6.79 1.15 1.92
## score3=2 23 12 4.25 14.15 19.09
##
## Chisq= 19.6 on 2 degrees of freedom, p= 0.00005
fit <- survfit(Surv(dfs, dfs.status) ~ score3, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.height = 0.2,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Disease-free survival by Model 3 score"
)
# By response to treatment
tab1(score3)
## score3 : Model 3 score
## Frequency Percent Cum. percent
## 0 24 33.3 33.3
## 1 25 34.7 68.1
## 2 23 31.9 100.0
## Total 72 100.0 100.0
tabpct(score3, respon3, graph=FALSE)
##
## Original table
## Response to treatment
## Model 3 score Complete Partial PD Total
## 0 21 2 1 24
## 1 16 6 3 25
## 2 11 7 5 23
## Total 48 15 9 72
##
## Row percent
## Response to treatment
## Model 3 score Complete Partial PD Total
## 0 21 2 1 24
## (87.5) (8.3) (4.2) (100)
## 1 16 6 3 25
## (64) (24) (12) (100)
## 2 11 7 5 23
## (47.8) (30.4) (21.7) (100)
##
## Column percent
## Response to treatment
## Model 3 score Complete % Partial % PD %
## 0 21 (43.8) 2 (13.3) 1 (11.1)
## 1 16 (33.3) 6 (40.0) 3 (33.3)
## 2 11 (22.9) 7 (46.7) 5 (55.6)
## Total 48 (100) 15 (100) 9 (100)
score3 <- factor(score3); pack()
surv.table(score3, 60)
## DFS Rate lower upper P-value OS Rate lower upper P-value
## Model 3 score "" "" "" "0" "" "" "" "0.038"
## 0 "95.2" "86.6" "100.0" "" "86.3" "72.8" "100.0" ""
## 1 "81.5" "66.6" "99.8" "" "81.6" "66.8" "99.7" ""
## 2 "33.1" "16.5" "66.3" "" "54.0" "33.2" "87.6" ""
## "" "" "" "" "" "" "" ""
# Overall Survival
# Figure 1 ####
survdiff(Surv(os, os.status) ~ stagegp)
## Call:
## survdiff(formula = Surv(os, os.status) ~ stagegp)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stagegp=Stage I 4 0 1.01 1.013 1.089
## stagegp=Stage II 20 3 4.56 0.532 0.768
## stagegp=Stage III 35 9 7.53 0.287 0.578
## stagegp=Stage IV 13 3 1.90 0.634 0.739
##
## Chisq= 2.5 on 3 degrees of freedom, p= 0.5
fit <- survfit(Surv(os, os.status) ~ stagegp, data = .data)
ggsurvplot(fit, data = .data,
# Change x-axis limit and scale
xlim=c(0,120), break.time.by=12,
# Add risk table
risk.table = TRUE,
tables.theme = theme_cleantable(),
# Color palette
palette = "lancet",
ggtheme = theme_bw(),
title="Overall survival by stage",
legend.labs=levels(stagegp)
)
