1. health and fitness
2. disease
3. cancer

Statistical Analysis of data from Breast
Cancer patients
MSc Dissertation
By Thomas Hunter Smith
School of Mathematical Sciences
University of Nottingham
I have read and understood the School and University Guidelines on plagiarism. I
confirm that this work is my own, apart from acknowledged references
Abstract
This report uses Survival Analysis techniques to analyse data on Breast Cancer patients
and attempts to determine factors which are associated with survival. Using the Nottingham Prognostic Index as a guideline and applying Kaplan-Meier, Log-Rank and Cox
Proportional Hazard Models; four main prognostic factors were identified to be directly
related to survival.
Contents
1 Introduction
1.1 Background . . . . . . . . . . . . . . . . . . . . . .
1.1.1 How did this particular project come about?
1.2 Methodology . . . . . . . . . . . . . . . . . . . . .
1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Breast Cancer and the NPI . . . . . . . . . . . . .
1.5 How it works..... . . . . . . . . . . . . . . . . . . . .
1.6 Geography of Breast Cancer . . . . . . . . . . . . .
1.7 The Evolution of Statistics in Medicine . . . . . . .
1.8 The Nottingham Prognostic Index . . . . . . . . . .
2 Data and Variables
2.1 Introduction to the Data . . . . . . . . .
2.1.1 A note on computing..... . . . . .
2.2 Missing Data . . . . . . . . . . . . . . .
2.2.1 Still some missing data though!! .
2.2.2 Complete Case Analysis . . . . .
2.2.3 Single Imputation . . . . . . . . .
2.2.4 Multiple Imputation . . . . . . .
2.2.5 Maximum Likelihood Estimation
2.2.6 Dealing with the missing data . .
2.3 Exploratory Analysis . . . . . . . . . . .
2.3.1 Relationships between Variables .
2.3.2 External Factors . . . . . . . . .
2.3.3 Method of Detection . . . . . . .
2.4 Summary of this Chapter . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
6
6
6
7
7
8
9
10
10
12
.
.
.
.
.
.
.
.
.
.
.
.
.
.
19
19
20
20
21
22
22
23
24
25
29
37
38
41
43
3 Survival Analysis
44
3.1 Survival Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 The Hazard Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3 Applications of Survival Analysis . . . . . . . . . . . . . . . . . . . . . . . 47
1
3.4
3.5
3.6
3.7
Censoring . . . . . . . . . . . . .
Purpose and Presumption . . . .
Non-parametric Survival Analysis
3.6.1 Kaplan-Meier Estimates .
3.6.2 Greenwoods formula . . .
3.6.3 Life Table Method . . . .
3.6.4 Nelson-Aalen Estimators .
Log-Rank Test . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
49
51
51
52
55
61
61
62
.
.
.
.
.
.
.
66
67
69
70
70
78
78
78
Research
. . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
82
82
83
85
4 Cox Proportional Hazard Models
4.0.1 Log Likelihood test . . . . . . . . .
4.0.2 Other Methods of Model Checking
4.1 Model Building using R . . . . . . . . . .
4.1.1 Model Building . . . . . . . . . . .
4.2 Index . . . . . . . . . . . . . . . . . . . . .
4.3 Model Checking . . . . . . . . . . . . . . .
4.3.1 Schoenfeld Residuals . . . . . . . .
5 Summary, Conclusion & Further
5.1 Summary and Conclusion . . .
5.2 Further Research . . . . . . . .
5.3 Acknowledgements . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Appendices
86
A Some further Interesting findings
87
2
List of Figures
1.1
Lymphatic System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12
2.13
2.14
2.15
2.16
2.17
Lymphatic System . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Estimated Survival Curves of Chemo-status amongst patients . . . . .
Scatter Plot of NPI against Survival . . . . . . . . . . . . . . . . . . . .
Scatter Plot of NPI against Survival with Fitted Survival Lines . . . . .
Scatter Plot of NPI for patients dying from Breast Cancer . . . . . . .
Scatter Plot of NPI for patients leaving the Study for Other Causes . .
Distribution of Survival for Patients Dying from Breast Cancer . . . . .
Distribution of Survival for Patients leaving the Study for Other Causes
Distribution of Survival for Patients still Alive at the end of the Study
Distribution of NPI scores across Centres (1-6) . . . . . . . . . . . . . .
Distribution of Tumour Size for Centre 1 . . . . . . . . . . . . . . . . .
Percentages of different Types of Tumour found in Centre 1 . . . . . .
Scatterplot of Tumour Size against Survival with fitted linear model . .
Histogram of age of Patients at 1st Operation . . . . . . . . . . . . . .
Histogram of age of Patients at 1st Operation . . . . . . . . . . . . . .
Distribution of Tumour Size for patients detected by screening . . . . .
Distribution of Tumour Size for patients not detected by screening . . .
3.1
3.2
3.3
Kaplan-Meier Curve computed from Table 3.6.1 . . . . . . . . . . . . . .
Kaplan-Meier Curves of Uncategorised NPI Scores . . . . . . . . . . . .
Kaplan-Meier Curves of NPI Groups with 5-year and 10-year Survival
Times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Kaplan-Meier Curves of Tumour Size (cm) categorised into 3 groups . . .
Kaplan-Meier Curves of different types of Tumour . . . . . . . . . . . . .
Kaplan-Meier Curves for different age groups . . . . . . . . . . . . . . . .
Kaplan-Meier Curves for Endocrine Therapy . . . . . . . . . . . . . . .
Kaplan-Meier Curves for Menopausal Status . . . . . . . . . . . . . . . .
3.4
3.5
3.6
3.7
3.8
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
26
28
30
32
33
34
34
35
35
36
37
38
39
40
41
42
42
. 54
. 56
.
.
.
.
.
.
57
58
59
60
62
64
4.1
4.2
4.3
Kaplan-Meier Curves of Histological Grades of Tumour, as well as of the
three factors involved in calculating the grade: Pleomorphism P, MitosisM
and Tubule Formation T . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Estimated Survival curves for 7 groups of new variable Total=T+P+M . . 72
Plots of scaled Schoenfeld residuals against transformed time for each covariate in the final CPH model . . . . . . . . . . . . . . . . . . . . . . . . . 80
A.1 Histogram of Survival for patients dying within 10 years of prognosis with
predicted survival of 85% . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
A.2 Kaplan-Meier Curves of NPI with 95% CI . . . . . . . . . . . . . . . . . . 89
A.3 Kaplan-Meier estimates of Estrogen Receptor Status . . . . . . . . . . . . 90
4
List of Tables
1.1
1.2
NPI Scores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Co-variates and their characteristics . . . . . . . . . . . . . . . . . . . . . . 18
2.1
2.2
2.3
Overview of the Breast Cancer Centres . . . . . . . . . . . . . . . . . . . . 21
Different Combinations of NPI . . . . . . . . . . . . . . . . . . . . . . . . . 29
NPI 5-Year Scores against actual Scores from Centre 1 . . . . . . . . . . . 30
3.1
3.2
Censoring Percentages across all centres . . . . . . . . . . . . . . . . . . . 51
Computation of Kaplan-Meier curve from section of Survival . . . . . . . . 53
.
5
Chapter 1
Introduction
”The Nottingham Prognostic Index (NPI) is an accepted prognostication tool in the management of breast Cancers. The latest overall 5-year and 10-year Breast Cancer survival
has been reported to be 85% and 77%”,
- BreastCancerResearch.com
1.1
Background
The aim of this report was to develop a thorough understanding of the branch of statistics
known as Survival Analysis and its application in Primary breast Cancer research. This
would be approached by evaluating the creation, function and application of the Nottingham Prognostic Index as well as applying methods used in Survival Analysis to a set of
data on Breast Cancer patients in an attempt to identify factors associated with survival.
1.1.1
How did this particular project come about?
Having not been particularly inspired by any of the dissertation topics kindly provided
by the University, I took it upon myself to seek an alternative subject. My passion
for medical research prompted me to look within the NHS. I managed to make contact
with Dr. Chan, a consultant at the Oncology Department at Nottingham City Hospital
and Professor Graham Ball, head Biostatistician at the John Van Geest Cancer research
centre, Nottingham. I was lucky enough to come to an arrangement which would give
me the chance to explore new statistical techniques on a real data set with information
on 16,965 Breast Cancer patients. I hope to take full advantage of this great opportunity
and deliver a final report which will include a comprehensive analysis of the data, and
perhaps include findings which may be of use in future Breast Cancer research.
6
1.2
Methodology
One of the first decisions made was to cut down the analysis to just one section of the
data which contained information on 2,203 patients. A comprehensive analysis of the
full 16,965 patients with the fluctuating inclusion and exclusion of certain co-variates and
data was unrealistic given the time frame of the report. Several different methods were
considered when dealing with the missing data, Little & Rubin (2002), including imputation and complete-case analysis - with the final consensus being that single imputation
is the most appropriate approach, implemented by an algorithm constructed in statistial
software package, R. The analysis of histograms was used to investigate distributions of
co-variates, and other visualisation techniques were used to examine the behaviour of the
data. Kaplan-Meier curves were used to compare survival times amongst the different
outcomes of a variable and the existence of significant differences between groups was
tested using log-rank. Proportional Hazard Models, Cox & Oakes (1982), were used to
identify which factors are significant in predicting a patients survival. Model checking
techniques included the Log-Liklihood ratio test, Kendall (1973), and the use of residuals,
Schoenfeld (1981) to identify violation of Proportional Hazards.
1.3
Results
Following the analysis, four prognostic factors remained significant in determining survival: Tumour Size, Histological Grade, Lymph Nodes (all of which make up the original
prognostic index) and Lymphovascular Invasion (LVI). Lymph Node involvement was
deemed to be the most important factor in determining prognosis for Breast Cancer. Despite prior thoughts to the contrary, this report found that a patients menopausal status
remains an insignificant factor, as was proven through both the log-rank test and Proportional Hazard models. The same was also true for age, although a re-categorised age
variable hinted there is evidence to suggest that these factors are not entirely insignificant
in determining prognosis. Time dependent factors such as recurrence and metastasis were
influential in determining survival, but were not considered while generating a Prognostic
Index using Proportional Hazard models, as they are unlikely to be present as prognostic
factors at the time of diagnosis. The same can be said of therapy variables. The presence
of a certain type of therapy on an individual is dictated by the biology of the tumour, thus
a patient’s survival chances must also be dictated through the biology of tumour. Rather
this than a patient’s prognosis being based on a set of information which is based on
another set of information, potentially leading to bias, loss of information and confusion.
It appears that the methods which are used to measure tumour size aren’t always fully
accurate. While measurements were being made to the nearest millimetre, an unusual
number of measurements were rounded to the nearest 0.5 of a centimetre. While this
7
rounding is for the most part compensated by the rescaling of tumour size in the Nottingham Prognostic Index, it did not fully remove the inevitability that on occasion there will
be rounding which will effect the prognostic index in a way that will alter the predicted
survival of a patient. Amongst the most significant findings was that a patient was up
to 20% more likely to survive from Breast Cancer if they have regular screening than if
they don’t. The 10 different types of tumour which existed, within the set of patients
analysed, did not prove to be an important factor in survival, although a patient with
Tumour Type 1 seemed to have smaller survival % than a patient with one of the other
nine types.
1.4
Breast Cancer and the NPI
”It is already the top Cancer in women both in the developed and the developing world,and
in many countries, it’s becoming more common.”,
World Health Organization, 2011
Breast Cancer makes up almost 1/5 of all Cancers found in women. In 2010, in the UK
alone, over 11,000 women lost their lives as a result of invasive Breast Cancer. As such, the
UK government continues to invest millions of pounds each year to aid and contribute to
various different areas of research, which share the common goal of prolonging a diagnosed
patient’s life for as long, and as pain free, as possible. Statistics is one of these areas.
• In 2010, 49,564 women and 397 men in the UK were diagnosed with invasive Breast
Cancer.
• In 2010 there were 11,633 deaths from invasive Breast Cancer in the UK.
• The most common symptom of Breast Cancer in men is a hard, painless lump that
develops on one of the breasts.
• In 2005-2009, 85 percent of women in England survived their invasive Breast Cancer
for five years or more.
• Breast Cancer is often thought of as only affecting women. Men can also develop
it, although it’s much less common. It affects just one in every 100,000 men in
England.
• In 2010, 5,765 women and 26 men in the UK were diagnosed with non-invasive
Breast Cancer.
8
1.5
How it works.....
Everybody is well aware of Cancer and the devastating effects it can have on people and
their lives, but surprisingly little people know HOW the Cancer can become fatal. In
tumours that grow inside of some vital organs, the very presence of the tumour itself
can interfere enough with the organ’s function that it eventually imparts lethal effects on
that organ. A good example of this would be a brain tumour. A brain tumour causes
tremendous amounts of localised swelling and it doesn’t have to be a large thing in order
to begin to crush vital brain tissue, leading to coma and death. Tumours of the lung can
plug off airways in a manner that leads to respiratory compromise, failure, and then death.
Cancers of the colon can obstruct the intestine or perforate, and both of these problems
can be fatal to the patient. Some tumours that arise from hormone secreting tissue such
as the adrenal glands or thyroid can actually cause such a great release of hormones that
it overwhelms the patient, leading to a fatal state. Tumour growth requires a tremendous
amount of bodily energy expenditure. This is fast growing tissue and, as a result, uses
up all of the normal energy stores and building blocks for healthy tissue such as protein.
Cancer risk increases with age:
Most Cancers are caused by DNA damage that accumulates over a person’s lifetime. Cancers that are directly caused by specific genetic faults inherited from a parent are rare. But
we all have subtle variations in our genes that may increase or reduce our risk of cancer
by a small amount - World Health Organisation
Breast Cancer, isn’t going to cause death due to the presence of the tumour in the breast
because the breast isn’t a vital organ. Breast Cancer can mestastasize to other organs
such as the liver and cause enough local destruction there to kill by way of organ damage, but even this is rare. Breast Cancer most often kills by spreading to many sites
in the body and creating enough tumour burden so that the patient is physiologically
overwhelmed, and they become cachectic and are susceptible to secondary illnesses such
as infections which ultimately lead to their demise. As the total mass of cancerous tissue
spreads throughout and the body accumulates multiple kilograms, it doesn’t matter if
it’s interfering with specific organ functions any more. The tumour burden itself starts
to make the patient ill. They begin to waste away. At that point, it’s a lot like having
AIDS or other such devastating diseases. The body becomes unable to combat even the
most minor illnesses it becomes exposed to. In these cases, patients succumb to secondary
processes like pneumonia or other infections. Not an experience that I would wish upon
anybody. So there are several different ways in which a Breast Cancer can develop and
harm the body.
9
1.6
Geography of Breast Cancer
There are surprisingly large differences in the incidence of Breast Cancer across different
parts of the world. Belgium has the highest rates at 109 out of 100,000 people (0.109%),
with other developed countries such as Australia with a 0.084% incidence rate. Compare
these to developing countries such as India (0.023%) and Mongolia which has as little
as 8 people per 100,000 diagnosed. No one has an answer to why this is the case one
train of thought is that women have children and begin breast feeding at an earlier age in
developing countries which is thought to reduce a womens chances of developing breast
cancer. But this still does not explain relatively large variations amongst countries where
the average age a women gives birth are very similar; Lifestyle, inheritance patterns, there
are many potential reasons. On the other hand, recovery rates are higher in developed
countries. ”Breast Cancer survival rates [range] from 80% or over in North America,
Sweden and Japan to around 60% in middle-income countries and below 40% in lowincome countries,” according to the WHO. The aim of this report is to draw upon the
medical data taken from 16,000 patients diagnosed with breast cancer and, through the
use of Survival Analysis and Proportional Hazard models, construct a model which will
accurately model a patients 5-year survival score. This is created with the intention
of providing the patient with as much accurate information as possible regarding their
Cancer, to help them decide upon the best course of action for themselves. There does,
of course, already exist a model which does exactly this.
1.7
The Evolution of Statistics in Medicine
It wasnt until around the 17th century that statistics was used in medicine. The idea of
statistical analysis itself was a relatively fresh and emerging concept. And indeed its use
in medicine would continuously breathe new life into medical research up right up until
the present day.
Prior to this very successful marriage, both concepts had lived very different lives.
Statistics had only really been around for a couple hundred years, while medicine itself
has been a primal instinct of humans ever since the dawn of civilisation. Prehistoric man
treated illness as a spiritual event, and would perform a whole manner of ritualistic methods to try and remove these demons. Popular techniques included Trephination (piercing
holes on the skull to ventilate these spirits) and the application of herbal medicines. For
some reason the former of these has fallen out of favour with most modern healthcare
systems, while the latter has proved more popular, and is still widely used in areas such
as Homeopathy. The evolution of medical practises was specific to the part of world it
orginated, with Europe, China, India and Native America all having their own remedies
to medical issues. It was not until the introduction of Medical Schools that medicine
really began to develop as a science, the first one of these is claimed to have been in
10
Salerno, Italy in the 13th Century. This allowed brilliant minds to work together, and
the renaissance led to revolutionary changes in the theory of medicine. In the fteenth
century, Vesalius repudiated Galens incorrect anatomical theories and Paracelsus advocated the use of chemicals instead of vegetable medicines. From the sixteenth-century
development of the microscope by Janssen and Galileo to the seventeenth century theory
of the circulation of blood (Harvey), scientists learned about the actual functioning of the
human body. The eighteenth century saw the development of modern medicines with the
isolation of foxglove (digitalis) by Withering, the use of inoculation (against smallpox) by
Jenner, and the postulation of the existence of vitamins (vitamin C, antiscorbutic factor)
by Lind.
Our knowledge (or indeed the invention and development) of Statistics on the other hand,
endured a far slower and less eventful growth spiral. Of course medical care has always
been an important part of our existence, but up until a few hundred years ago no-one
would have thought that mathematical theory could be used to help sustain and improve
human life. One of the birth places of statistics and probability was said to have been
born around 300AD when ancient greeks wanted to describe the outcome of a game of
dice between Zeus, Hades and Poseidon. Its theory developed through to the 13th century with the book of Abacus by Fibonacci and the 15th century with Luca Pacccioli
with his basic principles of probability. Cardano, Galileo, Pascal, Fermat, Bernoulli (and
his nephew) and Gataker all had significant inputs. These ideas were further refined by
scientists including Huygens (1657), Leibniz (1662), and Englishman John Graunt (1660)
wrote further on norms of statistics, including the relation of personal choice and judgment to statistical probability. Graunt in particular provided some of the first evidence
of Statistics being employed in medicine when he categorized the cause of death of the
London populace using statistical sampling. James Lind, a Royal Navy surgeon, carried
out the rst recorded clinical trial in 1747. In looking for a cure for scurvy, he fed sailors
aficted with scurvy six different treatments and determined that a factor in limes and oranges (subsequently found to be vitamin C) cured the disease while other foods did not.
His study was not blinded, but as a result (although not for 40 years)limes were stocked
on all ships of the Royal Navy, and scurvy among sailors (limeys) became a problem of the
past. By the 19th century medical research received a boom with the fast development
of technology to aid the understanding of the human anatomy. This gave birth to Case
studies which would use statistics to prove whether or not a treatment was beneficial.
Most research done before the twentieth century was more anecdotal than systematic,
consisting of descriptions of patients or pathological ndings. Many if not most of the
Statistical tests which we are familiar with today were established in the 20th century
by Sir Ronald Fisher, Archie Cochrane and Austin Bradford Hill, the first of whom is
considered by many to be the father of modern statistics. Epidemiological studies and
statistical modelling became increasingly popular, and lead to continued advancements
11
such as accelerated time failure and proportional hazard models. Statistics developed
for the purposes of medical research became commonplace across medicine, including of
course, Cancer research.
The publication by J. Mayer (2002) has a particularly fascinating take on the evolution
of medicine and statistics.
1.8
The Nottingham Prognostic Index
The Nottingham Prognostic Index (NPI) is used to determine prognosis following surgery
for Breast Cancer. Its value is calculated using three pathological criteria: the size of the
lesion; the number of involved lymph nodes; and the grade of the tumour. The formula
is shown in(1.1)
NPI = N + G + {0.2 × S}
(1.1)
Where:
• S is the size of the index lesion in centimetres (ie the Size of the Tumour)
• N is the number of lymph nodes involved: 0 = 1, 1 − 3 = 2, > 3 = 3
• G is the grade of tumour: Grade I =1, Grade II =2, Grade III =3
The resulting scores from this formula are then divided up into scores which would determine a patients predicted 5-year survival score. These scores are displayed in Table 1.8
These scores would be taken in to account when determining the patients prognosis, and
Score
5-Year Survival (%)
≥ 2.0 to ≤ 2.4
> 2.4 to ≤ 3.4
> 3.4 to ≤ 5.4
> 5.4
93
85
70
50
Table 1.1: NPI Scores
provide them with accurate information on their cancer, and their survival prospects.
The Nottingham Prognostic Index was originally derived from clinico-pathological correlation of a large series of symptomatic Breast Cancer patients. It is likely that it will need
to be revised in the near future to allow its effective application to screening patients.
- Breast Pathology on the Web
12
4
NPI
5
6
7
Figure 1.1: Lymphatic System
●
●●
●
●
●
●
●● ●
● ● ●
●
●●
●
●
●●●●● ●
●
●
●
●
●
● ● ●● ●
●
● ●
●
●●●
●●
●
●● ●●● ●
● ●
●● ●
●
●
●
●
●
●
●
●● ●
●
●
● ●
●
●
●
● ● ●●●
●
● ●
●
●
●
● ●●●
●●
● ●
●
●
●
●
●
●● ●●●
● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ●●
●●
● ● ●
●
●
● ●
●●
●
●●●
●●
●●●
●● ●
●
● ●●
●●●●
●
●●
●
●
●
● ● ●
●●
● ●●●●●●
●● ●●●●●
●
●
●●
● ●
●●● ●●
● ●
●
●●● ● ●
●●
●
● ● ●
●●
●
●
●
●
●
●
●
●
●
●● ● ●● ● ● ●●●
● ●
●
●●●●●
●●●●● ●
●
● ●●
●
● ●● ●
●
●
●
●
●
●●
●
●●●●
● ●● ●●● ●●●● ●
●
●● ● ●●
● ●
●
●
●
●
●
●
●
●●
●
●●
●●●
●
●●
●●● ●
● ●● ●●
●●●●●
●
●●
●
●
●
●
●
●●
●
●
●●
●●
●
●●●
●●
●
●●●●●
●● ●
●
●●● ● ● ● ● ●
● ● ●
●●
●●●●
●● ●●
●
●
●●
● ●● ●
●
●
●
● ●
●
●
●
●
●
3
●
●
●●●●
●●
●
●
●
● ●
●
●
2
●
0
●
●
● ●
●
●
● ●
●●
50
●
●
●●
●
●●
●
●
●
● ●●
●
●
●●●
●
●
● ● ●●
●
●
●
●
●
●
●
●●
●
●●● ●
●
● ● ●● ●
●
●
●
●
100
150
●● ● ●
200
Survival (months)
Originally formulated over 30 years ago, many Cancer centres across the world have used
the above 5-year survival score derived from the NPI as one of their primary sources to
determine prognosis following surgery. While this index has proved very successful over
the years in terms of pairing up a patient to the correct surgery, scientists and statisticians
are unanimous in their opinion that there is plenty of room for improvement. An example
of such an updated version would be the development of NPI+; a new updated version
of the old index. One fact which is widely accepted across all prognostic index’s is that
the most important factor in determining prognosis is the patients Lymph Node Status.
This status is based on the number of nodes from the lymphatic system involved in the
tumour (pictured) While this revised index is considered by most to be a more proficient
at prescribing the correct surgery, there is, as always, continued room for improvement.
While conducting research for this report, many online forums were encountered regarding
the correct prognostic index, with many sufferers received contrasting answers from different places. It seems that there are still large variations in the methods used to determine
prognosis out there, even with in the UK. Unfortunately, no unique and widely accepted
method has been established and this leads to a great deal of worry across women diagnosed with Breast Cancer. This uncertainty could be detrimental to patients chances of
13
survival due to the delays it will cause in their choice of treatment. It seems ridiculous,
but this amount of disparity draws comparisons to the inconsistencies of prehistoric medical practise mentioned earlier.
Hopefully one day, our improved knowledge of the biology of Breast Cancer and the advancement of ever more sophisticated statistical techniques will see a Unified Prognostic
Index developed which will optimise a patients survival time. The best way of understanding mathematical function is to see it in action, so let us have look at some example
of a patients Nottingham Prognostic Index taken from our study:
Patient ID 4231 had a size of 2.1cm, a histological grade of 3 and no Lymph nodes involved. The simple calculation of 0.2 × 2.1 + 1 + 3 yields a prognostic index of 4.42. From
the index, this provides a an estimated 5-year survival score of 70%. This would have provided the patient and the NHS with valuable information to aid them when deciding the
next best approach for treating the Cancer. As the 5-year survival score was relatively
middle ground in relation with other NPI scores, there might have been some difficult
choices with regards to the kind of therapy to go with. As it turned out, the patient
underwent Radio Therapy and was still alive at the age of 60, almost 9 and a half years
after the Prognostic Index was assigned. This is an example of the NPI working well;
the score would have suggested a relatively intensive treatment to maximise the patients
survival, and it paid off.
Patient ID 4612 had a tumour size of 4.8cm, 4 lymph nodes involved and a histological
grade of 3, yielding a very high prognostic index of 6.96. This comprehensively gave the
patient a 5-year survival score of 50 percent. This would mean that the patient would
require chemotherapy or run a high risk of not surviving more than 5 years. This index
score once again worked well, as the poor prognosis it administered prompted the patient
to go ahead with chemotherapy and was still alive 7 and a half years after the study began
at the age of 65.
Patient ID 2271 was looking in reasonably good shape at the beginning of the study and
was assigned a Prognostic Index score of 4.4, which is equivalent (according to the NPI)
to a 5-year survival score of 70 percent. The patient was told she stood a good chance
of living past 5 years, and could perhaps decide against undergoing chemotherapy. As
a result, the patient did not undergo any kind of chemotherapy treatment and unfortunately died at 58 only 3 years after the prognosis. This is an example of the Index going
wrong the patient may well have been misled with her prognosis and skipped therapy
which while being severe, may well have saved her life.
The final patient (ID 2758) provides an example that the Nottingham Prognostic Index
was used by the patient to make a concious decision about how to proceed. In the three
examples above patients were not forced into any kind of therapy assigned from the NPI,
14
the index merely provides them with information on their chances of survival and recommends them with an appropriate plan of action which will hopefully prologue their lives
for as long as possible. This particular patient was given a very poor prognosis with an
index score of 6.9, rating their 5-year survival at 50% or less. This NPI score would normally imply a patient undergoes Chemotherapy. But on this occasion data suggests the
patient weighed up the options and decided to stay away from the very undesirable effects
of chemotherapy and live out the rest of her live in as stress free a manner as possible.
She ended up dying 3 years later at the age of 41. These moral issues are very hard to
comprehend at the best of times, let alone attempting to explain them using mathematics.
So this report shall for the main part focus away from these delicate matters and focus
on delivering a prognostic index which will give the patient as much accurate information
as possible and make their choices based on well-grounded, truthful medical evidence.
Many more patients would need to be run through the NPI to develop a good idea of how
well the index works, but already its usefulness is evident. The application of the index,
and its effectiveness will be looked at in greater detail later on.
The process of constructing a prognostic index is an arduous one. First and foremost
the statistician would be looking for as large and complete a data set as possible, and of
course with the appropriate information. Information on the age of the patient at first
diagnosis and their current status (living or dead) is simply not enough. In the case of
this analysis, as much information on the biological and surgical status of each patient is
essential. The current Nottingham prognostic index considers only 3 factors associated
with Breast Cancer, but it has long been known that it can be improved, most likely with
the addition information from other areas. The dataset analysed in this study ticked most
of these boxes: It had information on more than 16,000 patients suffering from invasive
Breast Cancer. There did, inevitably exist sections of missing data, but nothing that
could not be managed. There was information on all biological factors suspected to relate
to breast cancer as well as information on all of the surgery and therapy the patients
went through. The information was taken over a set period of time, so the patients would
not be monitored until death but only to the end of the study. It is wholly unfeasible to
conduct this level of monitoring over a long period of time for more than a few patients.
Thankfully there are techniques available to deal with this. With the exception of a few
small cases of missing, there is as much information as one could hope for at this stage.
In addition to these are Time Dependent Variables: Distant Metastases, Local Recurrence, Contralateral Breast Cancer and Regional Recurrence. These were omitted from
the Model building process, the reason for which shall be explained later. With all of
this information at hand the analysis could begin. This would roughly follow the process
displayed below:
15
1. Conduct a thorough exploratory analysis of the data. Through the use of graphical
analysis and statistical tests, develop a thorough understanding of the co-variates
and their behaviour eg. using tests to determine significance and histograms to
understand distributions
2. Once the basic knowledge of the data has been established, begin to explore intereactions and correlations within the data set, which will provide useful information
for later on when constructing a statistical model. Pearson correlation,chi-squared
test
3. Use appropriate techniques to solve the problems regarding missing data, these techniques are context dependent but the main aim is to maintain as much of the data
as possible while minimizing the bias in the process
4. Begin applying survival analysis to the data, such as calculating and plotting KaplanMeier curves, to understand differences in survival amongst variables and their outcomes Begin to sub categorise variables using log ranks tests
5. Fit Cox Proportional Hazard models to the data and slowly but surely remove variables which have little to no impact on survival, and eventually have a model which
only contains variables (and any potential interactions and necessary transformations) which effect survival from breast cancer
6. Assess the validity and accuracy of the model by running new data through it, and
measure its expected outcomes against observed outcomes. Plot residuals and test
to ensure model does not not reject the assumption of proportional hazards Wald,
Likelihood and Score. Check diagnostics: for violation of assumption of proportional hazards, for influential data and for non-linearity in the relationship between
the log hazard and the co-variates
16
7. Using the coefficients from the model, construct a new prognostics index with new
5-year survival scores
In order to complete as effective an analysis as possible we must be able to weigh up all
potential factors effecting prognosis, and whenever possible, disregard any prior beliefs
which may distort potential truths which spill out from the study. It is often the case
(and indeed was the case at several points during the study) that a piece of information
pops up which completely contradicts previous inferences on the topic. More often than
not, this is a mistake on the statisticians part or perhaps a mis-collection of data. But
every so often this isn’t the case, and the findings might be very interesting and open new
windows for discovery. Data analysis should be entered with as open a mind as possible
(within reason), and should never disregard findings which lie in the face of those previous,
unless they can be categorically disproved. Employ tried and tested methods but also try
new approaches- History tells us that people have been blinded to the truth for centuries,
simply by not even questioning a long standing and respected statement. Breast Cancer
is one of the foremost areas of medical research, and is always at the tip of a biomedical
statisticians tongue. This report may not be that elusive cure but hopefully will add a
little piece to the jigsaw which scientists from all over the world are trying to complete.
17
Variable (R-code)
Description
Menopausal Status
meno
1=Pre-Menopause, 2=Post-Menopause
63% patients post menopausal
Lymph Nodes
lymphnode
1={0 Lymph Nodes Involved},
2 = {1 − 3}, 3 = {3 <}
Tumour Type
type
Nominal Categorical
10 different types
Age at first Operation
age
Range from 18-70
Mean age 54
Tumour Size(cm)
size
Continuous Variable
Range={0.1, ..., 4.9}, Mean=1.85cm
Method of Detection
detection
1=Screening 2=Other
Screening→ regular checkups
Lymphovascular Invasion
LVI
1=Cancer cells in Blood vessels
2=Not present in Blood Vessels
Tubule formation T
> 75% =1, 10 − 75% =2, < 10% =3
Mitosis M
Mitotic Count: 1=0-9,
2=10-19,3=20+
Nuclear Cells 1=Ok
2=Bad, 3=Worse
1 = {3, 4, 5}, 2 = {6, 7}
3 = {8, 9}
Indicator Variable
1=Therapy 2=No Therapy
”” ””
Pleomorphism P
Histological Grade (T+P+M)
grade
Chemotherapy
chemo
Radio Therapy
radio
Axillary Surgery
axillary
Endocrine Therapy
endo
”” ””
”” ””
Table 1.2: Co-variates and their characteristics
18
Chapter 2
Data and Variables
2.1
Introduction to the Data
The dataset provided by the NHS contains information on 16,965 women diagnosed with
Breast Cancer taken from 10 centres around the world. This data set provides a set of
information on each diagnosed patient, mainly the biological and surgical points of interest
which are considered to be relevant in this study. A brief overview of the variables and
what is going on is displayed in Table 1.8. A detailed description of these variables, in
the context of Breast Cancer is provided in the glossary towards the end of this report:
When dealing with such a large data set, it is very easy to become overwhelmed by the
sheer size of the task which lies ahead. It is very difficult to know where to start, which
way you should approach the analysis, what goals you should be aiming for, and so on.
- A thorough understanding of Survival Analysis and the Nottingham Prognostic Index
would be needed before beginning any serious analysis. Since R would be used extensively throughout the analysis, the anticipation and interpretation of its output would
be crucial. Analysis is often a trial and error game, but that is all for nothing if you
cannot tell between a meaningful result and a meaningless result
- If a dataset of this size was ever to be analysed with some useful findings under the given
time frame, some compromises would need to be made. Would a thorough analysis of
all 16000+ patients be possible? Would all of the variables be maintained ; some of
which surpassing 50 percent missingness?
- The data needed to be as complete as possible, so techniques would need to be identified to handle the missing data. Ideally these would maintain as much information as
possible without introducing any serious bias
- Efficiency in the choice of methodology, it is no use spending weeks using one technique,
then deciding that it is no good; I would need to be strategic with my selection
19
- In addition, there are never any guarantees than any real breakthroughs will come
about, so knowing when to stop is important-no research is ever pointless!
2.1.1
A note on computing.....
The analysis for this research project was conducted solely on Microsoft Excel and R,
both of which are excellent packages. In Excel, the built in package ”Data Analysis” was
used extensively for the exploratory analysis, and was very useful for summarising and
producing basic plots of the data. The R package ”survival” was used throughout the
Survival Analysis. This has many features, including inbuilt functions which can quickly
produce survival curves as well as fit and check Proportional Hazard Models.
2.2
Missing Data
It is never a massive surprise to come across missing data when dealing with large datasets
in any field; be it engineering or analytics or medicine, for whatever reason, the existence
of missing data is almost a forgone conclusion. Having considered the magnitude of the
task ahead, and the great volumes of analysis to be performed, it was decided early
on that the focus would be on just one centre. Given the time frame of this study,
a comprehensive study across all 10 centres, all of which with having vast differences
in the information provided (such as the variables described and missing data), would
realistically be infeasible. Looking at the varying information on the patients provided
across all 10 centres, it was important to choose a centre with minimal amounts of missing
data. This goes without saying, but perhaps more important consideration , would be to
choose a centre where the missing data (each one of the 10 centres had varying counts
of missing values, spread across all variables), it was insignificant enough to not distort
the analysis in any significant way. After careful consideration and preliminary analysis,
it was finally decided that the study would be focused on Centre (1). The main reasons
being as follows:
• The Number of patients was large enough for there not be any doubt over the
accuracy of the study (Which may have been an issue if, say, centre (2) was looked
at) and was digestible enough for R to handle the analysis, which again may have
been an issue if ,say, centre (9) was looked at.
• Using prior information on the variables which may play an important role in survival, centre 1 has an almost full data set, with almost complete information across
all variables considered important for the purposes of constructing an accurate and
effective CPH model
20
Centre
1
2
3
4
5
6
7
8
9
10
No. Patients
2203
1308
562
2391
1321
604
1539
1491
4136
1410
Mean NPI
4.07
3.23
3.83
3.45
3.66
3.92
3.63
3.83
3.45
4.12
Table 2.1: Overview of the Breast Cancer Centres
• The information the censoring of patients is almost complete, as a pose to most
of the other centres, where there are vast gaps in the information, which would of
course make matters a lot more complicated
2.2.1
Still some missing data though!!
One of the main reasons Centre (1) was chosen as the focal point of the study was the
relatively small ammounts of missing information it had when compared to other centres.
Survival analyses based on a data collection process which the researcher has little control
over are often plagued by problems of missing data. Deleting cases with any missing data
will result in information loss and usually results in bias, while many analytic procedures
that retain this information in some form underestimate the resulting uncertainty in
parameter estimates and other output. Despite this, there were still some gaps in the
information on many of the key variables:
• Several of the variables (eg. Parafin Block, Frozen tissue) which were monitored in
other centres were not monitored in centre (1). For the purposes of this study, these
variables shall be disregarded. In perhaps future studies, these variables can be
analysed across the other centres, and cross referenced with the results from centre
(1), to unearth their potential importance
• There some missing points in most of the remaining variables. Ideally the solution
to this problem would not involve completely removing patients with these missing
points from the study.
Several papers were covered in the search for the best approaches to handle missingness
in a survival analysis dataset. The correct/most effective approach is obviously context
dependent but Little and Rubin (2002) stood out for its clarity and general content.
21
These approaches are more applicable for higher levels of missingness, and only one of the
variables had higher than 5% missingness. So some amount of exploration into different
methods would be useful. Such an exploration would (i) undoubtedly be useful for future
data analysis, such as going on to analyse the other centres with greater chunks of information missing, (ii) provide further insight into the data itself; it is not unusual to make
interesting discoveries as bi-products of other analysis. There were 4 main approaches:
(a) Complete Case analysis (CC) (b) Single Imputation (SI) (c) Multiple Imputation (MI)
(d) Maximum Likelihood Estimation (MLE) Like almost all statistical techniques, these
have their benefits and drawbacks.
2.2.2
Complete Case Analysis
A direct approach to missing data is to exclude them. In the regression context, this
usually means complete-case analysis: excluding all units for which the outcome or any
of the inputs are missing. There are a few different ways to do this; list wise deletion (an
entire patient would be removed from the study if they had any missing data points) and
pairwise deletion ( the patient would be removed if there were data missing from any of
the variables used in the study). Immediately these strike as nave solutions for a study
with as many covariates as this. Generally speaking, there are two main problems arise
with complete-case analysis:
1. If the units with missing values dier systematically from the completely observed
cases, this could bias the complete-case analysis
2. If many variables are included in a model, there may be very few complete cases, so
that most of the data would be discarded for the sake of a simple analysis
The small amount of missing values are well spread across the data, and the application
of CC would see almost 11 % of patients being removed from the study, so this approach
has not been deemed viable.
2.2.3
Single Imputation
Albeit a relatively simple approach, single imputation is a very easy and efficient way to
replace missing values in a data set. Once again there are several paths one can go down
with this methodology:
1. Matching: This involves identifying a variable which correlates well with the one
with missing points, and comparing data points, so for example we know that
menopausal status has a strong relationship with age so a missing value in either
could be inferred from the other
22
2. Hot-deck imputation involves sorting a dataset according to any of a number of
variables, thus creating an ordered dataset. The technique then finds the first missing value and uses the cell value immediately prior to the data that are missing to
impute the missing value. Another method is known as cold-deck imputation, which
substitutes a value (with similar correlations) taken from another cold dataset of
the same nature
3. Rth nearest neighbour: a similar approach to matching and is essentially an algorithm for imputing missing data. Rather than comparing variables of single observations, it looks at numerous observations to conjure a value to impute. This is very
useful when looking at genetic arrays where there are 10s of thousands of missing
points
4. Mean substitution: substitutes the mean of the known data of the variable for the
missing data
A consideration for these methods is dealing with categorical data such as deciding which
value to take when that calculated lies in between the categories. For example when testing
mean substitution on the variable Did the patient receive Endocrine Therapy? the mean
value from the known data was 1.57 (1=Yes, 2=No) so how would we decide which answer
to give to the 52 patients with no information? Around 80% of the variables in the dataset
are indicator or categorical, so this is of obvious concern. In addition, mean imputation
ruins the correlations among variables because it changes the way two variables co-vary.
That is, if two variables move in tandem then using mean imputation to fill in the missing
values destroys the natural relationship between the predictor and outcome variables since
the imputed values will all be the mean of the predictor and remains static for all missing
values while the outcome still varies independently.
2.2.4
Multiple Imputation
Becoming increasingly popular with the continued development of computing power, multiple imputation methods are the most popular and generally most effective method of
dealing with missing data. In order to deal with the problem of increased noise due to imputation, Rubin and Little(1987) developed a method for averaging the outcomes across
multiple imputed data sets to account for this. The way this works is that imputation
processes similar to stochastic regression are run on the same data set multiple times and
the imputed data sets are saved for later analysis. Each imputed data set is analyzed
separately and the results are averaged except for the standard error term (SE). The SE
is constructed by the within variance of each data set as well as the variance between
imputed items on each data set. These two variances are added together and the square
23
root of them determines the SE, thus the noise due to imputation as well as the residual
variance are introduced to the regression model.
*statistical noise is a term that refers to the unexplained variation or randomness that is
found within a given data sample or formula
There are many, many different techniques when implementing multiple imputation, sometimes they can be tailored to suit a set of data exclusively and sometimes they can be
used across a whole range of datasets. Below are a few examples:
(i) multiple imputation without inclusion of the outcome (MI-)
(ii) multiple imputation with inclusion of the outcome (MI+)
(iii) Multiple Imputation using regression switching with predictive mean matching
(MICE)
(iv) multiple imputation using regression switching with predictive mean matching (MICE
PMM)
(v) multiple imputation using flexible additive imputation models
2.2.5
Maximum Likelihood Estimation
Maximum likelihood involves the estimation of parameters based on the observed data
using Bayesian estimation techniques. So MLE can use information already known about
a dataset, and use that to estimate a parameter of interest, say missing data. A likelihood function is the probability or probability density for the occurrence of a sample
configuration x0 , x1 , ..., xn given that the probability density f (x | θ) with parameter θ is
known.
L(θ | x) = f (x1 | θ)f (x2 | θ); ...f (xn | θ)
n
Y
=
f (xi | θ)
(2.1)
i=0
The Maximum Likelihood Estimator (MLE) θˆ can be obtained by maximising L(θ | x) as
a function of θ. As such it is usually far more convenient to work with the log-likelihood:
l(θ | x) = logL(θ | x)
n
1X
f (xi ; θ)
=
n i=1
(2.2)
∂
ˆ There are also more developed methods
l(θ; x) = 0 To obtain a value for θ.
Then set
∂θ
which use iterative procedures for parameter estimation, namely the E-M Algorithm.
24
The iterations alternate between performing an expectation (E) step, which creates a
function for the expectation of the log-likelihood evaluated using the current estimate for
the parameters, and a maximization (M) step, which computes parameters maximizing
the expected log-likelihood found on the E step. The ebook by Millar (2011) gives a more
modern take of MLE and how it has evolved along with computing.
Maximum likelihood estimation is equivalent to multiple imputation in most cases though
maximum likelihood is easier to implement than multiple imputation thus usually more
preferred, but multiple imputation has slightly more situations in which it is applicable.
2.2.6
Dealing with the missing data
Having already disregarded Complete Case analysis, the choice remained between single
and multiple imputation. When necessary, multiple imputation can be very effective but
it can also be very complicated. It was decided that in this instance it was not necessary,
as all but one of the variables are less than 5 % missingness. Any potential extra bias
which would be brought on by using single imputation will be minimal, and would not
effect the outcome and accuracy of the analysis.
In the Endocrine therapy example mentioned above, rather than replace the 52 missing
values with the value 2 as the average was 1.57, a function in R was used to randomly
insert the values based on the values already taken for this variable. So in this example the
random.impfunction randomly inserted 32Nos and 20 Yess into the missing data points.
random.imp <- function (a){
missing <- is.na(a)
n.missing <- sum(missing)
a.obs <- a[!missing]
imputed <- a
imputed[missing] <- sample (a.obs, n.missing, replace=TRUE)
return (imputed)
}
This algorithm worked extremely well for the majority of variables missing points, however there were a few which required more attention:
Internal Mammary Nodes
Of the 1564 patients accounted for in this variable, 98.15% 0 Internal Mammary Nodes,
with 1.66 % having 1 and 0.29 % having 2. With around 30 % missing, this was the variable had by far the most missing information (with all others being less than 5 percent
25
Figure 2.1: Lymphatic System
700
120.00%
600
100.00%
500
# of Patients
80.00%
400
60.00%
300
40.00%
200
20.00%
100
0
0.00%
3
4
5
6
7
8
9
More
T+P+M
missingness), imputing new values would be a relatively arduous task. With this being the
only variable with substantial missing information, conditions are ideal for non response
weighting to be applied. This would be a long and relatively fruitless task, considering
the fact that it correlates strongly with number of Axillary nodes positive. Due to this
high level of correlation, the decision was made to
Tubile Formation, Pleomorphism and Mitosis
These three factors are used to calculate the final Histological Grade of the tumour.
Twenty five patients had all three of these missing at once, and a further 5 had just Mitosis
missing. Fortunately, every single patient had their final histological grade, so it wasn’t
difficult to impute values such that they didn’t violate the criteria. For example, one
patient had information on the final grade, as well as Tubile formation and Pleomorphism
both being graded 3. The histological grade was 2, so the only value Mitosis could
take would be 1, otherwise there would be inconsistency in the model. The sub grade
information for T, M, P for those patients with all three missing was imputed under the
assumption that:
• those with a histological grade of three would automatically have a grade of 3 for
Tubile formation and Pleomorphism
26
• those with a histological grade of one would automatically have a grade of 1 for
Mitosis and Pleomorphism
Having manually filled in the points to match this criterion, the random.imp function
could be used to randomly insert the rest of the missing data for T,P and M.
Recurrence
The problem with these recurrence variables is that their potentially hazardous effects
don’t usually become present until some time after the patient has had initial surgery.
This means recurrence is a time-dependent covariate, which has the potential to make
the model building process far more complex. It is clear from the data that the recurrence variables have a large say in the patient’s survival, so it is important that they are
taken into account when modelling survival. Unfortunately the fact that they remain
anonymous until sometime until after the patient’s first prognosis means that they will
essentially be useless in a Prognostic Index as it is unlikely that any sign of them will
appear until some months and years after the patients first 5-year survival score is given.
There does exist methods which accommodate Time dependent variables into a model,
but in this situation they are not required: while it will be possible how the existence
of recurrence effects survival, it is unlikely that a model will exist which can accurately
predict if recurrence will occur. As such the time dependent variables which involve recurrence were removed from the study.
Chemotherapy and the 9s....
An interesting problem arose when it was unveiled that 59 (2.7%) patients had observation
which contained 9s in place of where the standard 1’s and 0s would be for Chemotherapy.
Figure 2.2 shows that this third outcome has even lower survival times than those who
underwent chemotherapy. One possible answer to this is that patients who had this outcome did undergo chemotherapy, but ceased the treatment as their cancer had developed
to a point that survival seemed unlikely. The choice was made to keep these data points as
they are, and treat chemotherapy as an ordinal categorical variable. This decision makes
sense in that there is a clear difference in the survival times for these three outcomes.
Later research using the log-rank test confirmed a clear separation in survival between
the three outcomes.
Summing up this section, there was a substantial amount of lists and a,b,c-ing which
explains the great variety of methods which exist for handling missing data. As was
mentioned earlier, the right approach is always context dependent but, even in matters
regarding a single issue there are often several approaches to take. Another major con27
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 2.2: Estimated Survival Curves of Chemo-status amongst patients
0.0
Chemotherapy
No Chemotherapy
Unknown Therapy
0
50
100
150
t
28
200
sideration are the time restraints put in place. In the case of this report the choice of
approach was not complicated by this factor due to the low levels of missingness which
required only the simple matter of single imputation. The application of methods of multiple (regression based )imputation would probably have resulted in reduced biased, but
give the low levels of missing data and the given time restraints, single imputation was
deemed adequate.
2.3
Exploratory Analysis
A long time was spent exploring variables, their relationships with each other, distributions and effect on survival. Considering the behaviour of the original NPI against
survival provided a reasonably good starting point. We can see that the prognostic index
scores are segregated into five clear groups, which makes good sense if we break down the
combinations of the score in Table (2.3). S = 0.2×Tumour Size ≈ 0.36 ± 0.12, based on
the median values for tumour size, is assumed be consistent for each combination.
This is reflected in Figure 2.3, with the patients clustered around the NPI regions of
Group
(1)
(2)
(2)
(3)
(3)
(3)
(4)
(4)
(5)
Lymph Node Status
1
1
2
2
3
1
3
2
3
Histological Grade
NPI Score
1
≈ 2.36 ± 0.12
2
≈ 3.36 ± 0.12
1
”” ””
2
≈ 4.36 ± 0.12
1
3
”” ””
2
≈ 5.36 ± 0.12
3
”” ””
3
≈ 6.36 ± 0.12
”” ””
Table 2.2: Different Combinations of NPI
(2.24,2.48), (3.24,3.48), (4.24,4.48),(5.24.,5.48),(6.24,6.48) and is even more evident in
the distribution of the NPI scores in Figure2.10. The size of these regions is down to
the first and third Quartiles of Size multiplied by the 0.2 as noted in the index. The
data points lying outside of these regions were patients with either exceptionally small or
exceptionally large tumour sizes.
The majority of the observations lie in the groups (2) and (3) of the index. What may
be of interest now, and potentially useful for looking at constructing a new index later
on is to see how these five groups differ in survival. For groups (1) and (2), the majority
29
4
2
3
NPI
5
6
7
Figure 2.3: Scatter Plot of NPI against Survival
●●● ●
●
●●
●●
●
●
●●
● ● ●
●
●
●●
●
●●●
●
●
●●
●● ●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●●●●● ●
●
●
●● ●
●
●
● ●
●
●
●
●
●
●● ●●● ●●●
●
●
●
●● ●
●
●●
● ● ●●
●●● ●
●
●
●
●
●
●
●
●
● ●● ●●
●
●
●●
●●
●
●
●
● ● ●●● ● ●●
● ●●
●● ●
●●
● ●●● ● ●
●●●
●●
● ●●●●
●●
●
● ●
●●●
●●
●
●● ●●
●●
●●
●●●●●
● ● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●● ●●● ●
●●
●● ●●
● ●
●●●
●●
●●
●●
●●
●●●
●
●
●
●● ● ●●
●●●
● ●●● ●
●● ●
●
●
●
●
●●●
●
●
●●
●
●●
●
●
●
●
●●
●●
●
●
●
●●
●
● ●
●●●
●●●
●
●●●● ●●
●
●
●
●●
●● ●
● ●●
●●●
●●
●●●●
●●●
●
●
●
●●●
●
●
●●
●
●
●
● ● ●●
●
●●●
●
●
●
●●
●
●
● ●●●
●●
●●
●●
●
●●
●●
●●●
●
●
●
●●●●
● ●
●●●
●●●●
●●
●●● ●●●●●●●
●●
●
●
●
●●●● ● ●
●
●
●
● ● ● ●● ● ●● ●●
● ●●
●●●
● ●
●●
●●●● ●● ●● ●● ●●
● ●●
●
●
●
●
●
●
●●
●●
●
●
●
●●●●
●
●● ●
●● ●
●
● ●
●
●
●
●
●
●● ●●
●● ●● ●●● ● ●● ●●
● ● ●
●●
●●
●
●
●● ●
●
●
●● ●
●
● ● ●●
●
●●
●
●●
●●●●
●
●
●
●● ●●
●
●●
●●
●
●
●●
●
●
●●
●●
●●
●●
●●
●
●●
●
● ●● ●
●
●
●
●
●●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●●●
●
●
● ● ●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●●●
●
●●●
●●
●
●●
●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●●● ●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●● ●●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●
●● ●●●
● ●
●
● ● ●●
●
●●●●
●●
●
●●
●
●
●
●
●●●●
●●
●●●● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ● ●●
● ●
●
● ●
●●●● ●
●
●
●
●
●
●
● ●● ●
●
●
● ● ● ● ●● ●●
● ●●
●
●●
●
●●
●●
● ●●
●
●
● ●●●●●
●●
●●
●●
●
●
●
●
● ●
●
●
●
●●
●●●●●
●
●
●
●
●
●●●
● ●●●
●●
●
●
●● ●
●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●●
●
●●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●●
●
●
●●
●
●
●
●
●
●●
●●
●●
●●
●
●●
●
●
●
●
●
●
●
●●●
●●● ● ●●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●●
●
●
●●
●●●●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●●
●● ● ● ●
●
● ●●● ●
●
●
●
● ●
●
●
● ● ●● ●
●
●●
●
●
●●
●●● ●
●●
●
●
●
●●
● ●
●
● ●
●
●●
●
● ●● ●●●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●
●
●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●●
●●
●
●
●
●
●
●
● ●●●●
●●
●
●●
●●
●
●●
●●
●
● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●●●
●●
● ●● ●●
● ●● ●
●●
●
● ● ●●●
●
●
●●
●
●● ●
●
●
0
50
100
●
150
200
Survival (months)
of the patients had a survival time of more than 80 months, and while the same applies
for (3), there is an clear increase of patients surviving < 80 months. A measure of the
accuracy and performance of the NPI against the 2203 patients from Centre (1) can be
easily obtained. Table (2.3) shows the 4 different categories of NPI scores, the number of
patients and the actual survival percentages of each group in Centre (1).
It is immediately transparent that there are huge differences between the index score
5-Year Survival (NPI Score)
90% (≥ 2− ≤ 2.4)
85% (> 2.4− ≤ 3.4)
70% (> 3.4− ≤ 5.4)
50% (≥ 5.4)
Number of Patients
% Survival Actual
321
96.4% (307 patients)
468
95.1% (445 patients)
1107
86.2% (954 patients)
308
61% (188 patients)
Table 2.3: NPI 5-Year Scores against actual Scores from Centre 1
and the actual score - is this intentional? Looking at the different categories of 5-year
survival, these patients do not match very well with the predicted survival percentages. The lower scores (< 2.4)for the NPI matches well with insignificant difference
30
(pairedt − test, α = 0.05) however the three subsequent categories are less encouraging.
For the next two stages, the prognostic predicted survival differs significantly(α = 0.05)
with 10% and 11% difference respectively, with more patients surviving beyond the 5-year
mark than the index anticipated. And the final category anticipated 10% more people
surviving then actually did. Additionally 85% survive longer than 5 years, which emphasises how the vast majority of patients have low range NPI scores. It is initially surprising
to see that less than 50% of the patients survive for less than 10 years, but this is likely to
be a figure severely distorted by censoring, meaning that the survival times for this data
take into account patients who were still alive at the end of the study, from which point
their age cannot be monitored. Thus it that age is used for the analysis. Since many of
these studies will last less than 10 years, the true percentage survival at 10 years is tricky
to establish. The 5-year survival is far easier to clarify, as almost all studies of this nature
will last at least 5 years. So we have seen there are significant inconsistencies, and you
would expect something as developed as the NPI to be more accurate there are several
possible reasons for this.
Were possibly not enough patients tested on? There seem to be different levels of survival
across centres Centre (1) seems to have reasonably high levels
Perhaps the index is designed to be under predict survival by around 10% to make
patients think more cautiously about their choices. I would be interesting to hear the
truth, and if indeed this is the truth, is it morally right for NHS to distort information
for the patients in such a way?
Or perhaps the NPI is that inaccurate and in desperate need to revising?
Interesting thought......when counting up the patients surviving less than 5 years after diagnosis there was quite a high intensity of patients dying slightly after the 5-year
survival time, especially 6-7 years as can be seen in Figure A.1. This shows the distribution of time of death for patients with NPI score > 2.4to3.4 ≤ dying before 10 years.
This is one of a few examples where 5 years as time point of reference for survival seems
inappropriate. On this evidence, were this particular prognostic index to be used, either
the survival % for each group should be increased by 10% or the time point be increased
to 7 years. If this was to be performed to the scores from centre 1, the actual scores would
be insignificantly different (α = 0.05) to the NPI scores. Figure 2.4 demonstrates how the
predefined categories segment the data. The vertical red line represents the 5 year mark,
and the blue horizontal lines mark the different categories which determine survival. So
as we are moving upwards through the sections we are expecting the number of points on
the right to decrease and on the left to increase (in propotional terms for that particular
category). While this is the case, it does not fit in with the levels of survival stated by
the NPI. The green lines represent the actual times which would need to be assigned
to this dataset if the prognostic index was to maintain the survival percentages for its
31
4
3
2
NPI
5
6
7
Figure 2.4: Scatter Plot of NPI against Survival with Fitted Survival Lines
●●● ●
●
●●
●●
●
●
●●
●
● ● ●
●
●
●●
●
●
●
●
●
●
● ●●● ●● ●●
● ●●
●
●
●●
●
●
●
●
●
●
● ●
●● ●●
●
●
● ●
●● ●●●
●
●
● ●●
●
●
●● ●●● ●●●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●●
●
●
●
● ● ●●● ● ●●
● ●●
●● ●
●●
● ●●● ● ●
●●●
●●
● ●●●● ●
●
● ●
●●●
●●
●●
●
●● ●●
●●
●
●●●
●●
●
●● ● ●●
●
●
●● ● ● ●●
●
●
●●
●●●● ● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●●
●●
●
● ●●●●●
●
●●●
●
●
●● ●
●● ● ●●
●●●
● ●
●
●
●●
●●●
●
●●
●
●●●
●
●
●
●
●●
●
●
●● ●
●●
●
●
●●
●
●
●●
●
● ●
●●●
●●●
●
●●●● ●●
●●●●
●
●●
● ●●
●●
●● ●
● ●●
●
●
●●●
●
●
●●
●
●
●
● ● ●●
●
●●
●
●●●
●
●
●
●●
●
●
●●
●●
●
●
●●
●●●
●
●
●
●●●●
● ●
●●●
●●●●
●●
●●● ●●●●●●●
●●
●●
●● ●●
●
●
●
●●●● ● ●
●●
●
●
●●
●
● ● ● ●● ● ●● ●●
●●
● ●
●●
●●●● ●● ●● ●● ●●
● ●●
●
●
● ●
● ●●
●●
●●
●
●
●
●
●
●● ●
●
●
● ●
●
●
●●●
●
●
●● ●●
●● ●● ●●● ● ●● ●●
● ●
●●
●●
●
●
●
●● ●
●
●
●● ●
●
● ● ●●
●●
●
●●
●●●●
●
●
● ●●● ●●
● ● ●●
●●●
●
●
●
●●
●
●
●●
●●
●●
●●
●●
●
●●
●
● ●● ●
●
●
●
●
●
●
●●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●●
●
●●
●●
●●●●●●
●●●
●●●
●
●
●
●●●
●●
●●
●
●
●
●
●
●
●
●
●●
●
●●
●●
●
●
●●
●
●●●●
●●
●
●●
●
●
●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●●
●
●●
●
●
●
●
●●
●
●
●●
●
●
●● ●
●
●
●
●
●●
●●
●●●
●
●
●●
●
●
●
●
●
●● ●●●
●
●
●
●●
●
●●●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●●● ●●●
● ●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●● ●
●●● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●●
●
●
●● ●●●●●●●
●
●●
●
●
●● ●
●● ● ●
●
●
●
●
●●●
●
● ● ●●●
●
●
●
● ●●
● ●
●
●
●
●
●●●
●
● ●
●
●
●
●
● ● ●● ●●
● ●● ●
●
●
●
● ●●
●
●
●
●●
●●
● ●●
●
●
● ●●●●●
●●
●●
●
●●
●●
●
● ●
●●●
●●
●
●●
●●●●●
●
●
●
●
●
●●●
●
●● ●
●●
●
●
●● ●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●●
●●●
●
●
●●
●
●● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●●●
●
●
●●●
●
●
●
●●
●
●●
●●
●●
●●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●●
●
●●
●
● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●●
●●●
●●
●
●
●
●
●●
●
●●●
●●●●
●●●
●●
●
●
● ●● ●
●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●●
● ●●● ● ●
●● ● ● ●●
●
●
●
●
●
●
●
●
●● ● ●
●
●●
●
●●
● ● ● ● ● ●●
●
●
●
●
●●
●
●
●
● ●
●
●●
●
● ●● ●●●●
●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ●●●
●
●
●●●
●
●
●
●● ●●
● ●●● ●
● ●●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●●
●
●
●
●●
●
●
●●●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ●
●●
●
●●
●
●
●
●●
●●
●●
●
●●●
●●●
●
●● ●
●
●
●
●
●●
●●●
●●●
●●
●
●●
● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●
●
● ● ● ●●
●
● ● ●● ●
●
0
50
●
100
survival
32
150
200
4
NPI
5
6
7
Figure 2.5: Scatter Plot of NPI for patients dying from Breast Cancer
●
●●
●
●
●
●
●● ●
● ● ●
●
●●
●
●
●●●●● ●
●
●
●
●
●
● ● ●● ●
●
● ●
●
●●●
●●
●
●● ●●● ●
● ●
●● ●
●
●
●
●
●
●
●
●● ●
●
●
● ●
●
●
●
● ● ●●●
●
● ●
●
●
●
● ●●●
●●
● ●
●
●
●
●
●
●● ●●●
● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
● ●●
●●
● ● ●
●
●
● ●
●●
●
●●●
●●
●●●
●● ●
●
● ●●
●●●●
●
●●
●
●
●
● ● ●
●●
● ●●●●●●
●● ●●●●●
●
●
●●
● ●
●●● ●●
● ●
●
●●● ● ●
●●
●
● ● ●
●●
●
●
●
●
●
●
●
●
●
●● ● ●● ● ● ●●●
● ●
●
●●●●●
●●●●● ●
●
● ●●
●
● ●● ●
●
●
●
●
●
●●
●
●●●●
● ●● ●●● ●●●● ●
●
●● ● ●●
● ●
●
●
●
●
●
●
●
●●
●
●●
●●●
●
●●
●●● ●
● ●● ●●
●●●●●
●
●●
●
●
●
●
●
●●
●
●
●●
●●
●
●●●
●●
●
●●●●●
●● ●
●
●●● ● ● ● ● ●
● ● ●
●●
●●●●
●● ●●
●
●
●●
● ●● ●
●
●
●
● ●
●
●
●
●
●
3
●
●
●●●●
●●
●
●
●
● ●
●
●
2
●
0
●
●
● ●
●
●
● ●
●●
50
●
●
●●
●
●●
●
●
●
● ●●
●
●
●●●
●
●
● ● ●●
●
●
●
●
●
●
●
●●
●
●●● ●
●
● ● ●● ●
●
●
●
●
100
150
●● ● ●
200
Survival (months)
set categories. From this visual demonstration we see that we would need to squeeze or
retract the blue lines enough so that the green line matches the red line. This would
require altering the margins to something more or less than 1 (or 2), and obviously. The
nature of this obviously requires on the distribution of the points. Were they uniformly
distributed, it would be easier to see how the index would needed to be altered to fit the
data.
In the case of centre 1, the data are gathered in 5 groups, but there are only 4 categories
for these to fall into. There is no immediate explanation for this, but it is a reminder that
the model building is just part of the process, and the subsequent index requires some
amount of visualising.
Figures 2.5 and 2.6 represent the difference in NPI scores for patients dying from breast
cancer those leaving the study from other causes. These can help visualize the possibility
of informative censoring.
After looking at the age distributions for Death from Cancer in Figure (2.7), leaving the
study other causes (Figure 2.8) and still alive at the end of the study (Figure 2.9), they
all have strong suggestion that they belong to a family of normal distributions. But by
33
Figure 2.6: Scatter Plot of NPI for patients leaving the Study for Other Causes
●
● ●
●
●
6
●
●
●
● ●●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
5
●
●
●
●
● ●
●
●
●
●
4
NPI
●
●●●● ●
●
● ● ●●
●
●
● ●●
●● ● ● ●
● ●
●
● ●
●
● ●
●
●
●
●● ● ●
●
3
●
●
●
● ●●
●
● ●● ●
● ● ●
●
●
●
●
●
●
●
●
● ● ●
●● ●● ●
● ●
●
●
2
●
●
0
●
●
●
●
●
●
●
●●● ●
● ●●
●
●
●
●
●
●●
● ● ●● ●
● ●
● ●
●● ● ● ●●
50
●
●
●
●
●
●
100
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●●
●●
●
●
150
●
200
Survival (months)
Figure 2.7: Distribution of Survival for Patients Dying from Breast Cancer
485 patients dying of from breast cancer
40
35
Frequency
30
25
20
15
10
5
0
Age
34
Figure 2.8: Distribution of Survival for Patients leaving the Study for Other Causes
173 patients dying from other causes
25
Frequency
20
15
10
5
0
40
44
48
52
56
60
64
68
72
76
80
84 More
Age
Figure 2.9: Distribution of Survival for Patients still Alive at the end of the Study
1545 patients still alive at the end of the study
140
120
Frequency
100
80
60
40
20
0
Age
35
Figure 2.10: Distribution of NPI scores across Centres (1-6)
Distribution of Nottigham Prognostic Index Scores across
patients from Centres 1-6
800
700
Frequency
600
500
400
300
200
100
0.1
0.3
0.5
0.7
0.9
1.1
1.3
1.5
1.7
1.9
2.1
2.3
2.5
2.7
2.9
3.1
3.3
3.5
3.7
3.9
4.1
4.3
4.5
4.7
4.9
5.1
5.3
5.5
5.7
5.9
6.1
6.3
6.5
6.7
6.9
More
0
NPI
comparing these histograms, information can be gathered which may give more insight
into the nature of the censoring. The similar patterns all three groups follow follow suggest that the reasons for those who left the study for other causes, may be linked to the
breast cancer. While these deaths might not come directly from breast cancer, they may
be directly associated with it.
Spreading across other centres, Figure (2.10) looks at the NPI of over 8,000 patients. It
can be seen that the patterns of the histograms follow a similar shape to the other groups
looked at. But Immediately it is evident that something is not quite right as the graph
is displaying a number of patients with NPI scores of less than the base rate of 2. This
will be the result of either missing information on the Tumour Grade or the Number of
Lymph Nodes involved.
Looking at Figure (2.11), once again a normal distribution comes to mind. On closer
inspection it becomes evident that there is a disproportionate frequency of patients with
a Tumours Sizes with factor 0.5, implying that for more than a few patients, this measurement is being rounded up/down to the nearest factor of 0.5, which is surprising. The
histogram suggests that it is possible to accurately measure the size of a tumour to the
nearest mm, yet it is obvious that less accuracy has been used in many instances. Given
that the size of the tumour plays a relatively large role in determining prognosis, hence
the choice of surgery, hence the chances of a patients survival, perhaps more time and
36
Figure 2.11: Distribution of Tumour Size for Centre 1
Distribution of Tumour Size amongst patients in Centre 1
250
Frequency
200
150
100
50
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
3.2
3.4
3.6
3.8
4
4.2
4.4
4.6
4.8
5
0
Tumour Size (cm)
effort could be invested in measuring the tumour size as accurately as possible. Having
said that, the potential dangers of this measurement strategy will be limited substantially
be the fact that Tumour Size is only 20% present in the final NPI. It is quite possible that
those involved in the measurement process were well aware of this fact. Were Tumour Size
to carry more power within a Prognostic Index, perhaps there would be more emphasis on
the accuracy. Additionally it seems more accuracy is used for measuring smaller tumours
(< 1.5cm), and as the size increases, so does the amount of rounding.
There were 9 different types of tumour found in the patients. Figure (2.12) shows that
the majority are of type 1, with type 11d also making up a large share of the patients.
With such a range of tumour types, instinct dictates that differentiation between types
may influence survival in some way. A pie chart is no way to measure this but later on,
differences in survival between types will be investigated using Kaplan-Meier estimators.
2.3.1
Relationships between Variables
With the range of types of variable in the dataset, there does not exist a single generalised method to detect relationships. Since there exists continuous data (eg Survival),
nominal categorical (eg Tumour type), ordinal categorical (eg Histological Grade) as well
as binary/indicator variables (eg Endocrine therapy; a great deal of care must be taken
in choosing the correct method for each situation. A rigorous procedure of comparing
37
Figure 2.12: Percentages of different Types of Tumour found in Centre 1
% Spread of Tumour Types
T12
1.54%
Missing
0.50%
T11d
24.33%
T8
0.22%
T1
58.24%
T7
5.44%
T3
6.72%
T6
0.18%
T5
1.50%
T4
1.32%
variables is not a necessity but in an ideal world a good working knowledge of the main
relationships would be attained before throwing the variables into a statistical model.
Having already had a look at the behaviour of the NPI against survival, it might be of
interest to explore its relationship against another continuous variable, say size, as shown
in Figure 2.13
While not being immediately striking, the data points are mainly populated around about
90-170 months in survival with Tumour size (< 2.5cm. And there is a trend appearing
after about 90 months of survival, the smaller the tumour size the longer the survival as
might expect. This is echoed by the line, which is the linear relationship between the two
variables. While this is heading in the right direction, it is certainly not a true relationship
between the variables as it ignores several key attributes to the data such as censoring.
2.3.2
External Factors
There is always the possibility that there are factors which are not included in the study
which could affect a patients survival. Censoring (which shall be looked at in greater
depth later on) takes into account situations where the patients survival has been compromised by factors unrelated to Breast Cancer. But there are times where the collection
of information for variables maybe inconsistent amongst patients. For example the time
38
5
Figure 2.13: Scatterplot of Tumour Size against Survival with fitted linear model
●
●
●
●
●
●
●
● ●
●
● ●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●●
●
● ●● ●
●●● ●
●
●●
●
●
●
●
●
●
● ●●●
●
●
●
●
●
●
●
●
● ●●
● ●
● ●
●● ●
●
● ●
● ●●●
●● ●
● ● ●● ● ●●●●●
● ●●
●
●
●
● ●●●
●
●●●●●
●
● ●
● ●●
● ● ●
● ●●
●
●●
●
●● ● ●●●
●● ● ● ●
●●
● ●● ●
●
●
●●
●
●
●
●
●
●● ●
●
●●●● ● ●● ●
●
●●●
● ●
●
●
●
● ●
●●
●●●
●●
●● ●
●●●●
●
● ●●● ●●● ●●
●
●●
●●
●●
●●
●●●
●● ●●●●
● ●
●●
●●
● ●
●●
● ● ●●
●●●
●●
●
●●
●
● ●●● ●● ●●●●
● ● ●●
●
●●●●●●
●●
●● ●
●
●
● ● ●
●● ●●● ●
●●
●●●
● ●● ●
● ● ●
●●
●
●
●
● ● ●●
●
● ●●
●
● ●●
●●
● ●
●
●
●
●
● ●● ●●
●●
● ●●●
● ●
● ● ●●
● ●●
● ●
●
●
●●●
● ●● ●●●●●
●●
●●
● ●●●
● ●● ●●●●
●
●●●● ●● ●●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●●
● ●
●●
● ●●● ●●
●
●
●
●
● ● ● ●●●
●
●
●
●
●●●● ●
● ●●● ●●
●
●
●●
●●●
●
●
●●●●
●●● ●●●●●●●
●
●●
●●●
●●●
●●
●
●
●
●●●●●● ●● ●
●●● ●●●
●
● ●● ●
● ●●
●
●●
●●
●●● ●
●●
●●●
●●
●
●
●●● ●●
●●●
●
●
● ●
●●
●●●●
●●●●
●
●●●●
●
●
● ●●
●●● ●
●
●
●
●●
● ●●
●
● ●●● ●
●
●
● ●
● ●
● ●● ●
● ●●
●
●● ●●●
●
●●●● ●
●● ●●
●●●●
●●
●●
●
● ●●●
●
●●
●●
●●
●●
●
●
●●
●
●
●●●
●
●
●
●●
●
●
●●
●
●●●
●
●●
●●
●
●●
●
●●
●
●
●
●
●
●●●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●●
●
●
●●
●
●●●
●● ●
●
● ●● ●● ●● ●
●●●●●●
●●●
●
● ●
●●● ●●●
● ●●● ●
●
●●
●● ●
● ●
●
●
●●●● ●●●●●● ●●● ●
●
●
● ● ●
● ● ●
●● ● ●●●
●●●● ●
● ● ●● ●●
●●
●
●
●●●
●●●●
●
●●
●
●●●
●●
●●
●
●●
●●
●
●
●
●●●●
●
●●
●
●●●●●
●
●●
●●
●
● ●●● ●●
●
●●●●● ●● ●
●
● ●● ●●●●●
●
●● ●●
●●
●
●
●●
●
●●
●
●
●
●●
●●●●
●
●
●●●●●●
●●
●●● ●
●●●
●
●
●●●
● ●●●
●
●
●
●
●
●
●● ●
●
●
●●●
●
●●●
●●● ● ●● ●●
●
●
● ●●
●
●
●
●
●●
●●
● ●●
●●●
●
●●●●
●●●●
●
●
●● ●●● ●●
●
● ●
●
●●
●
●●
●●
●
●●●
●● ● ●
●
●●●●●
●●●
●●
●●
●●●●
●●
● ●●
●●
● ● ●●
●●
●
●
●●●
●●
●
●●●
●
●
●
●●●
●
●●
●
●
●
●●●
●
●●●
●
●●●
●
●
●
●
●●
●●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●● ●
●
●
●
●
●
●
●●
●●
●
●●
●
●●
●●● ●●
●
●●
● ●
●●●
●●
●
●●●
●
●
●●
●●●●●
●
●
●
●
●●● ●● ●
● ● ●●●●●●
●
● ● ● ●●
● ●
●●
● ● ●● ●
●
●
●●
● ●●
● ●● ●●
●
●
●
●●●●
●
●
●
●
●
●●●
●
●
●
●
●●●
●●
●
●
●●●●
●
●
●●●
●
●●
●●
●●●●●
●
●●●
●
●●
●
●
●●●●
●
● ●● ●●●
● ●●
●●
● ●
●●●●
●●●
●
●●●●
●●
●●
●●●
●
●
●
●●
●●
●
●
●
●
●
●●
●●●
●●●●
●
●●●
●
●● ● ●●●
●● ●
●
●
●● ● ● ●● ●●
●●
● ● ● ● ● ●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●●
● ●●●
● ●●
●
●● ●
●
●
●●
●●●
●
●●
●●
●● ●●
●●●
●●
● ●●
●●
●
●
●●●
●
●●●●●
●
●
●●
●●
●●
●
●
●●
●●●●●●
●●
●
●
● ●●
●
●
● ●●●●
● ●● ●●
●●●
●
●
●
●
●
●
●
●●●
●●
●
●
●
●●
●●
●
●
●●
●
●
●●
●
●
●●●● ●●● ●● ●
●●●●● ● ●
●●●
●●
●●●
●
● ●●
●
●●
●
●●
●
●
●
●
●●●
●●●
● ●●
●●
●
●●
● ● ●●●●●● ●●●
● ●
●
●●● ●
●
●
●
● ●● ● ● ●●
●
●
●●●
●
●
● ●
●
●
●●
●
●● ● ●●● ● ●●●
●●●● ●
●
●
●
●●
● ●● ●● ● ● ●●
●
●
●●●●
● ● ●● ●● ● ● ● ●
● ●● ●
●●●
●
●
● ● ● ● ● ●
●
●
● ● ● ●● ●
● ● ●●● ●
●●●
●
●
● ●●● ● ● ●
●
● ●● ●
●
●
●●
●
● ● ●
● ●●● ●
●
●●
●●
●●
●●
● ●
●
●
● ●
● ● ●●
●
●
● ●●
●
●
●
●
0
1
2
Size(cm)
3
4
●
●●●
●
0
50
100
Survival
39
150
200
Figure 2.14: Histogram of age of Patients at 1st Operation
Frequency
Patients age at 1st Operation
200
180
160
140
120
100
80
60
40
20
0
Age of Patient at 1st Operation
between a person having cycles of chemotherapy may be fluctuate between patients due
to different waiting times in hospitals. Or the time taken between initial diagnosis and
the first operation may be affected by a factor not looked at in this study, such as delays
in readily available equipment/medication or the mental state of a patient. One of the
variables that may be an important factor is the date of diagnosis; however information
on this is unavailable. The closest information relevant to this is the date of the 1st operation, the distribution of which is displayed in Figure 2.14. While there might be time
variations between the date of diagnosis and the date of first operation for each patient,
it should be noted that the cancers progression will be dictated by the characteristics of
the tumour, all of which are taken at the same point for each patient. Figure 2.14.
One variable which provides insight into the stage of cancer for each patient going into
the study is Method of Detection This refers to the circumstances in which the tumour
was first detected and has two possible outcomes: Detected through Screening or Not
detected through screening. The patients who had detected a tumour through screening
would have, on average, a less progressed breast cancer. On the other those who were
diagnosed out with screening are likely to have further progressed cancer. This makes
intuitive sense, as the regular frequency of screening would help ensure that the tumour
was identified as early as possible. Patients not having regular screening check ups run the
risk of a tumour sitting unnoticed for a long period of time, bringing with it far greater
40
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 2.15: Histogram of age of Patients at 1st Operation
0.0
Screening
Other Method
0
50
100
150
200
t
risk. A way of measuring the effect that late diagnosis can have on a patient is to plot a
Kaplan-Meier curve showing survival. Figure 2.15 shows the different percentage survival
rates of patients over time for patients who’s tumour was detected by screening against
those who was not. This striking difference in survival for these two methods of detection
may suggest that it may well be a factor in determining survival, but then again cancers
severity is dictated by its biological characteristics, which represent themselves regardless
of the time they were detected.
2.3.3
Method of Detection
Exploring the difference in methods against other variables provides interesting results.
Take Figures (2.16) and (2.17). Comparing the histograms on tumour size of Screening
against not screening roughly show two normal distributions, however the tumour size
for screening is centred around 1.2cm, while the histogram for the patients who did not
undergo screening was centred around 1.8cm. So if a patients tumour was first caught out
in screening, they can expect to have a tumour 60mm less in size than if their tumour was
detected another way. It is certainly true that there will be a relationship between the
therapy a patient undergoes and their subsequent prognosis. This does not necessarily
41
Figure 2.16: Distribution of Tumour Size for patients detected by screening
Tumour Size Histogram - Method One (722 Patients)
70
60
Frequency
50
40
30
20
10
5
4.8
4.6
4.4
4
4.2
3.8
3.6
3.4
3
3.2
2.8
2.6
2.4
2
2.2
1.8
1.6
1.4
1
1.2
0.8
0.6
0.4
0
0.2
0
Tumour Size (cm)
Figure 2.17: Distribution of Tumour Size for patients not detected by screening
Tumour Size Histogram - Method Two (1481 Patients)
140
120
80
60
40
20
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
3.2
3.4
3.6
3.8
4
4.2
4.4
4.6
4.8
5
Frequency
100
Tumour Size (cm)
42
depend on the success of the therapy, but on the nature of the therapy. For example, a
patient undergoing Chemotherapy treatment is more likely to have a poor prognosis than
a patient undergoing hormonal treatment. These are two different processes likely to be
prescribed at different developmental stages of the cancer. Chemotherapy is considered
the most severe treatment, and is only used if all other treatments are not successful. So
a patient undergoing chemotherapy is likely to have worse survival chances than patients
undergoing other treatments such as radiotherapy and Endocrine therapy.
2.4
Summary of this Chapter
Exploratory analysis is a fairly systematic exercise, and the results displayed in this chapter are amongst the most interesting from many, many comparisons, investigations and
plots. With the frequent use of Scatter plots and Histograms; visualisation is often the
way to unravel any secrets lurking in the data. While the exploratory Analysis is more
about trial, error and instinct, handlign the missing data is a far more subtle exercise. In
addition to Little & Rubin (2002), Allison,P (2001) contains a well refined list of methods.
Now the missing data has been dealt with, and any hidden secrets within the variables
has been unearthed, now is as good a time as any to move on with the Survival Analysis.
43
Chapter 3
Survival Analysis
The term ”survival analysis” is a statistical approach designed to take into account the
amount of time an experimental unit contributes to a study. That is, it is the study of
time between entry into observation and a subsequent event. Less formally, a unit may be
a particular patient under study or a structure under observation while the event would
be the focal point of the study such as death. Survival Analysis made an immediate impact in medical research in particular. When observing units as unpredictable as humans
over a certain time frame, problems would arise such as patients leaving before the timeframe of the study had elapsed. This could be for a number of reasons, such as quitting
the study, moving too far away to continue monitoring or dying from reasons unrelated
to the topic under investigation. As a result of these unpredictable and uncontrollable
events, the researcher was forced to withdraw the units and all associated data from the
study. This would have been extremely frustrating as units which would have provided
valuable information to the study, have to be completely disregarded. As a result these
studies could never be entirely regarded as trust worthy, as there was simply too much
information missing.
Thankfully techniques called censoring enabled researchers to analyse incomplete data
due to delayed entry or withdrawal from the study. Survival analysis meant that studies
could start without all experimental units enrolled at start of observation time and could
end before all experimental units had experienced an event. This allowed each and every
experimental unit to contribute all of the information possible to the study considering
the amount of time the researcher was able to observe the unit. In a similar fashion to
those involved in the construction of the Nottingham Prognostic Index, this project shall
focus on the use of Survival Analysis.
44
3.1
Survival Function
Survival Analysis relies on the Survival Function S(t). This function is a property of any
random variable that maps a set of events, usually associated with mortality or failure of
some system, onto time. It captures the probability that the system will survive beyond
a specified time. Let T be a non-negative random variable representing the waiting time
until the occurrence of an event. For simplicity this report shall adopt the terminology of
survival analysis, referring to the event of interest as ’death’ and to the waiting time as
’survival’ time, but the techniques to be studied have much wider applicability. Assuming
for that T is a continuous random variable with probability density function f (t) and
cumulative distribution function F (t) = P r{T ≤ t} giving the probability that the event
has occurred by duration t.
<t+4t)
f (t) = lim4t→0 P (t≤T4t
(3.1)
The Survival Function is often referred to as the Complementary cumulative distribution
function. This can be seen below, with the Survival Function conventionally denoted S(t),
for some time t and T (usually representing time of death or failure) be a continuous
random variable with cumulative distribution function F (t) on the interval [0, ∞). Its
survival function is
S(t) = P (T > t)
Z ∞
f (u)du
=
0
= 1 − F (t)
(3.2)
Simply put, (3.2) shows the probability that the time of death is later than some specified
time t. This has the following properties:
• Every survival function S(t) is monotonically decreasing.The survival function must
be non-increasing: S(u) ≤ S(t)if u ≥ t. This property follows directly from F (t) =
1 − S(t) being the integral of a non-negative function. This reflects the notion
that survival to a later age is only possible if all younger ages are attained. Given
this property, the lifetime distribution function and event density (F and f below)
are well-defined; the survival function is usually assumed to approach zero as age
increases without bound, i.e., S(t) → 0 as t → inf, although the limit could be
greater than zero if eternal life is possible.
• The survivor function is smooth and right continuous
• The time, t = 0, represents some origin, typically the beginning of a study or
the start of operation of some system. S(0) is commonly unity but can be less to
represent the probability that the system fails immediately upon operation. In other
45
words one would usually assume S(0) = 1, although it could be less than 1 if there
is the possibility of immediate death or failure.
This function would be used for two main reasons. Firstly determining a patients probability of surviving to time t and secondly determining the % which survive survive to
time t.
The roots of survival analysis can be traced back to mortality tables from centuries ago,
however it was not until the second world war that the Survival Analysis of the modern
day began rapid development. Ironically, his new era was stimulated by interest in reliability (or failure time) of military equipment. It seems that it was the divine incentive
to kill rather that to save which really pushed forward this breakthrough. At the end of
the war these newly developed statistical methods emerging from strict mortality data
research to failure time research, quickly spread through private industry as customers
became more demanding of safer, more reliable products. As the uses of survival analysis
grew, parametric models gave way to non-parametric and semi parametric approaches
(all of which I shall discuss in greater detail later on) for their appeal in dealing with the
ever-growing field of clinical trials in medical research.
3.2
The Hazard Function
The Hazard function h(t) is, in short, used to measure the failure rate for very small
intervals (or continuous) time. This provides a crucial link between Survival analysis
and Proportional Hazard models, and assesses the instantaneous risk of demise at time t,
conditional on survival to that time:
h(t) = lim4t→0 P rd{t ≤ T ≤ t + 4t} | T ≥ te/4t
(3.3)
In words: the probability that if a unit survives to t, that unit will succumb to the event
in the next instant.In terms of the survival function:
f (t)
s(t)
d
= − lnS(t)
dt
= slope of − lnS(t)
h(t) =
= how many (in %) die at time t
(3.4)
This can be interpreted by considering an individual alive at time t. The chances of dying
in a small interval [t, t + 4t) are then given by:
qt =
S(t) − S(t + 4t)
≈ h(t)4t
S(t)
46
(3.5)
Which can be alternatively written as:
h(t) =
N o of deaths in [t, t + 4t)
N o person-years at risk in[t, t + t)
(3.6)
It should be noted that, unlike the survival function S(t), the hazard function is NOT a
probability so can take values of less than 0 and more than 1. Also, unlike the survival
function, the hazard function can instantaneously decrease OR increase (OR both) over
time.The modelling of survival data usually employs a hazard function or the log-hazard.
For example, assuming a constant hazard, h(t) = v, implies an exponential distribution
of survival times, with density function p(t) = ve−vt . In essence the Survival increases
with t as the Hazard decreases with t and vice versa.
3.3
Applications of Survival Analysis
Despite the assumptions that one may make from its name, Survival Analysis is not only
restricted to modelling cases concerning the actual transition of an organism from living
to death. It has a much broader application, using essentially the same methodology
across various different fields. Some examples of this include event-history analysis used
in sociology, reliability theory used in Engineering and duration analysis in economics. It
is not all that difficult to imagine how this topic is so easily transferable across various,
and in many ways very different areas. . So a question posed in Medical research such as
What proportion of a group suffering from an undesirable condition will survive past a
certain point? could be answered with exactly the same statistical techniques as a question
asked by an incompetent Engineer What proportion of these houses will collapse before
a certain point? Or a psychiatrist looking for an answer to the question What effect does
electroshock therapy have on a patients mental stability? would need only go as far as the
answer to an economists question What effect does David Cameron have on the United
Kingdoms economic stability? or perhaps even a combination of the two What effect does
electro-shock therapy have on David Cameron?
Perhaps these examples are a little far fetched, but never the less give you a taste of the
diversity and transferability of Survival Analysis. As such, it is used increasingly in areas
such as Political science, where these approaches can be extremely informative.
– In international relations; war duration, alliance longevity, peace duration, length
of trade agreements, etc
– In generalized politics; Government durations, length of civil wars, , confirmation
durations, political careers, bill co-sponsorship, policy innovation/adoption coalition
durability, etc
47
– In America, where the political structure is a lot more complex than simply Democratics Vs. Republicans; Campaign contributions (which, on the evidence of last
years US elections , may well be worth Billions of dollars), Congressional overrides
of Supreme Court decisions, etc
Or in other problems such as economics, sociology, anthropology,
– Length of a career (within a single role, including redundancies and promotions etc),
Strike durations
– Issues dealing with criminals and their employability etc, success in life based on
marital/childhood circumstances, in music, the chances of a band being able to
make a record label a sizeable profit
– Various health and development related issues
More generally, survival models are useful for any issue in which the phenomenon of
interest is a duration, and/or the response is the occurrence of a discrete event in time.
As with any statistical approach, process, method or model, there are certain rules which
apply, and with these rules often come problematic issues. But first lets have a look at
the characteristics of Time-To-Event Data
1) Take place over time
2) Discrete events (i.e., not continuous)
3) the unit may not experience the event (i.e., possibility of censoring)
Criteria for Survival Analysis
Unfortunately these characteristics bring with them several problems and boxes which
must be ticked if they are to function in the correct manner. Bringing the concept of time
into any scenario makes things a whole lot trickier as is the case with Survival Analysis;
- all of the duration data is strictly non-negative (unfortunately time travel has still eluded
even the sharpest of scientists)
- Also, the data are conditional: to survive to some time t, one must necessarily have
survived up to t-1 as well
- Additionally, observations are regularly encountered which have not failed yet (i.e.,
censored data)
Thankfully there does exist well exercised methods of dealing with these problems, which
shall be covered in greater detail in the next section.
48
3.4
Censoring
Censoring is a form of missing data problem which is common in survival analysis. A
thorough understanding of censoring and its role are crucial to completing a comprehensive
analysis of the survival data, and hence fitting the proportional hazard models.
• Left censoring a data point is below a certain value but it is unknown by how much
• Interval censoring a data point is somewhere on an interval between two values
• Right censoring a data point is above a certain value but it is unknown by how
much
• Type I censoring occurs if an experiment has a set number of subjects or items
and stops the experiment at a predetermined time, at which point any subjects
remaining are right-censored
• Type II censoring occurs if an experiment has a set number of subjects or items and
stops the experiment when a predetermined number are observed to have failed; the
remaining subjects are then right-censored
• Random (or non-informative) censoring is when each subject has a censoring time
that is statistically independent of their failure time. The observed value is the
minimum of the censoring and failure times; subjects whose failure time is greater
than their censoring time are right-censored
Ideally, in a medical scenario (involving patient monitoring based on life and death) such
as this, both the birth and death dates of a subject are known, in which case the lifetime
is known. If it is known only that the date of death is after some date, this is called right
censoring. Right censoring will occur for those subjects whose birth date is known but
who are still alive when they are lost to follow-up or when the study ends. If a subject’s
lifetime is known to be less than a certain duration, the lifetime is said to be left-censored.
It may also happen that subjects with a lifetime less than some threshold may not be
observed at all: this is called truncation. Note that truncation is different from left
censoring, since for a left censored datum, we know the subject exists, but for a truncated
dataset, we may be completely unaware of the subject. Truncation is also common. In a
so-called delayed entry study, subjects are not observed at all until they have reached a
certain age.
Left censoring of data can occur when a person’s survival time becomes incomplete on
the left side of the follow-up period for the person. As an example, we may follow up a
patient for breast cancer from the time of them being diagnosed. We may never know the
exact time that said patient will have had the tumour for.
49
Looking at the censoring in the context of our dataset, the ”birth” of our patient is now
considered the point at which they were first diagnosed with Breast Cancer. At the end
of the study, there are three considered and categorised outcomes which each patient shall
face: They shall either
- Have survived beyond the period of study (Often the most likely scenario, as realistically
studies cannot extend to decades, and the expected survival time of any diagnosed
patient far surpasses this)
- Have died as a direct or indirect result of having breast cancer
- Have left before the end of the study
The latter of these outcomes can be from a variety of reasons; the patient may have
decided (through poor health or otherwise) that they no longer wished to be monitored,
they might have (as a direct or indirect result of breast cancer) suffered other serious
illnesses which became more serious etc. The table 3.1 provides more of an insight into
the nature of the censoring in this particular study, where
0 = The patient was still alive at the final point of study
1 = The patient died as a result of Breast Cancer
2 = The patient has left the study for other reasons
For obvious reasons, 3.1 provides further justification of the decision to focus the study
on Centre (1). The nature of this categorisation may seem relatively straightforward, however its application and legitimacy divides opinions of statisticians around the world. The
main point of discussion being how to draw the line between ”related” and ”unrelated”
deaths? How can one be certain that the ”unrelated” deaths are not as an indirect result
of having a tumour ?
For example, following surgery a patient may become more susceptible to other potentially
fatal infections. Lets say a patient has undergone chemotherapy, which can significantly
repress a bodies immune system, and they later are diagnosed and killed with pneumonia.
For the case of the data monitoring and collection, this death would usually be classified as
death by other causes. Naturally these questions open up a whole new range of problems,
and unfortunately cannot be fully addressed in this report (due to limited information on
the patients, time constraints etc). As a result, this study shall neglect this eventuality,
and concentrate on analysing the censored data in the form it is presented.
While it is likely that left censoring is present within this dataset, and that some of the
patients may have had had the tumour for some ammount of time before they were diagnosed, for the purposes of this project, all patients will be considered to be right censored.
50
Centre
1
1-5
6-10
Patients
2203
7785
9180
0
70.1%
73.2%
72.8%
1
22%
12.8%
11.4%
2
7.9%
3.8%
3.4%
Missing info on Censoring
0%
10.2%
12.4%
Table 3.1: Censoring Percentages across all centres
The vast majority cases in Centre (1), and indeed the entire dataset are Right Censored.
3.5
Purpose and Presumption
The 3 main goals of Survival Analysis are:
1. To estimate and interpret survival characteristics, for example looking a Kaplan
Meier plots of our variables to see how they behave against survival
2. Compare survival for different groups, for example using the log-rank test to see
how different outcomes of a variable effect survival
3. To asses relationship of Explanatory variables to survival time using a Cox Proportional Hazards model
These all run under two main assumptions:
1. Independent Observations
2. Independent Censoring
Independent observations, speaks for itself, the outcome of one patient will in no way
affect the outcome of another. Independent censoring: Consider two random individuals
both alive at time t, such that the first individual is not censored and the second individual
is censored at some point before t. Censoring is called independent if these individuals
have equal chances of survival.
3.6
Non-parametric Survival Analysis
Once the time to event data has been collected, the first task is to describe it usually
this is done graphically using a survival curve. Visualisation plays a large role in Survival
Analysis and allows the appreciation of temporal patterns in the data. It also helps us to
identify an appropriate distributional form for the data. If the data are consistent with a
parametric distribution, then parameters can be derived to efficiently describe the survival
pattern and statistical inference can be based on the chosen distribution. Non-parametric
51
methods are used when no theoretical distribution adequately ts the data. There are
three non-parametric methods for describing time to event data:
• the Kaplan-Meier method
• the Life table method
• the Nelson-Aalen method
All three of these methods have their place in modern day Survival Analysis, although
Kaplan-Meier is generally accepted to be the most effective technique.
3.6.1
Kaplan-Meier Estimates
Kaplan-Meier estimators shall play a crucial role in determining the varying effects which
the variables have on survival. The Kaplan-Meier curves which these estimators yield
show a plot of percentage survival against time. They are often a used in survival analysis
provide a very digestible graphical output, which can be easily interpreted by those less
familiar with the statistical concepts at work. The Kaplan-Meier method is based on
individual survival times and assumes that censoring is independent of survival time (that
is, the reason an observation is censored is unrelated to the cause of failure). Kaplan-Meier
estimator of survival at time t is shown in (3.7).
S(ti )
=
Q
0≤j≤n
P r(T > ti | T ≥ ti )
= S(tj−1 )P r(T > tj | T ≥ tj )
nj − mj
for 0 ≤ t ≤ T
= S(tj−1 )
nj
(3.7)
Here tj , j = 1, 2, ....n is the total set of failure times recorded, T is the failure time with
distribution F and density f , mj is the number of failures at time tj , and nj is the number
of individuals at risk at time tj . Note that:
1. for each time period the number of individuals present at the start of the period
is adjusted according to the number of individuals censored and the number of
individuals who experienced the event of interest in the previous time period, and
2. for ties between failures and censored observations, the failures are assumed to occur
first
A note on Kaplan and Meier....... Both were students of the famous John Tukey.
In 1952, Paul Meier started working on the duration of cancer while at Johns Hopkins
University, Chicago, USA. Edward Kaplan later started working on the lifetime of vacuum tubes in the repeaters of sub-oceanic telephone cables while at Bell labs. They
independently submitted their research on survival times to the Journal of the American
52
Statistical Association, whose editor encouraged them to submit a joint paper, which they
did: Kaplan & Meier (1958).
Here is an example of survival times (in months) taken directly from our dataset, the +s
indicate censoring has occurred:
88 66 179+ 188+ 121 136 179 177+
So the uncensored survival times show that the patient has deceased from breast cancer
that long after the study began. The censored times indicate that the patient was still
alive and at that age once the study period had expired OR that they had died from
other causes during the same period. The table 3.6.1 below shows how a Kaplan-Meier
would be computed from the above survival times by taking into account censoring. The
resulting curve is shown in Figure 3.1 - notice the red line’s corresponding to the values
of S(tj ).
• t(j) = Ordered failure times
• mj = N o of patients at tj )
• nj = # failure times at t(j)
• qj = # censored in [tj , tj+1 )
• S(tj ) =
#survivingpastt(j)
n(j)
t(j)
mj
nj
qj
S(tj )
nj − mj
t(0) = 0 m(0) = 8 n(0) = 0 q(0) = 0 S(t(j) ) =
S(tj−1 )
nj
8−7
× 1 = 0.875
66
8
1
0
8
88
7
1
0
0.75
121
6
1
0
0.625
136
5
1
0
0.5
177
4
0
1
0.5
179
3
1
1
0.335
188
3
0
1
0.335
Table 3.2: Computation of Kaplan-Meier curve from section of Survival
53
0.0
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.1: Kaplan-Meier Curve computed from Table 3.6.1
0
50
100
t
54
150
3.6.2
Greenwoods formula
The KaplanMeier estimator is a statistic, and several estimators are used to approximate
its variance. Once such approximation is:
X
qj
ˆ j )) = S(tj )2
(3.8)
V ˆar(S(t
n (nj − qj )
t <t j
j
This result in (3.2) is known as Greenwoods formula. Its derivation may be questioned
because it conditions on the nj which are random variables, but the result is in the
spirit of likelihood theory, conditioning on all observed quantities, and has been justied
ˆ
rigourously. Peterson (1977) has shownq
that the K-M estimator S(t)
is consistent, and
ˆ − S(t)) converges in law to a Gaussian
Breslow and Crowley (1974) show that n(S(t)
process with expectation 0 and a variance-covariance function that may be approximated
using Greenwoods formula.
As already been seen, Kaplan-Meier curves are very good at giving an immediate idea of
the behaviour of different groups with survival, especially when the variable of interest
can only take 2-4 values. One of the good things about K-M curves are that they are
very quick and easy to interpret, and as such one can skim through large amounts of
variables/outcomes and inference behaviour in a small space of time.
Since NPI is a continuous variable and can take many values, there are many survival
curves on the Plot 3.2, which is quite frankly a mess. Not much information at all can be
taken from this so a way of organising the NPI data must be found such that it will not
result in the Survival curves looking like a box of string. If the NPI is intuitively broken
down into groups, then making sense of its relationship with survival is a possibility.
Spitting a continuous variable up into a categorical variable can be a delicate issue as
there are several points which need to be addressed:
- How many categories will the variable be group up into? Are there any trends in the
variable which can provide clues on the best way to do this?
- Shall the categories be of even interval or shall they be adjusted to cater for specific
behaviour in its continuous form?
- Are there methods in place to help with these decisions?
In the case of Figure 3.2, the variables were separated into 4 groups with the same intervals
as the index (In Table 1.8). These intervals are not of even length, so the scientists
behind the NPI decided that there were patterns in the Prognostic Indexs relationship
with survival which suggested that intervals of equal length were not appropriate. Figure
55
0.0
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.2: Kaplan-Meier Curves of Uncategorised NPI Scores
0
50
100
150
t
56
200
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.3: Kaplan-Meier Curves of NPI Groups with 5-year and 10-year Survival Times
0.0
5−Year Survival
10−year Survival
0
50
100
150
200
t
3.3 and shows the survival curves of the 4 groups of the NPI.
There are methods of checking whether these categories are significantly different from
one another but from just looking at the plots we can deduce information about the NPI
scores from this set of 2203 patients. The curves are in the right ordering as we might
expect, with 2 < S ≤ 2.4 being the group with highest survival and S > 5.4 being at the
bottom, where S is the a patients NPI score at the beginning of the study. There is a
small amount of overlap between 2 < S ≤ 2.4 and 2 < S ≤ 3.4 categories of the index,
but the different curves seem to reflect the chosen intervals of index well. If you refer to
A.2 in the appendix, there is minimal overlap of CIs once again with exception to the top
two. The 5-year survival times from this are displayed in Table 2.3. This shows that the
NPI was miles off with its 5-year survival times. Coincidently, it seems like the actual
10-year survival time is closer to the 5-year survival times postulated by the index.
Another continuous variable which might have more predictive power when split into
categories and is also currently useless in a KM-plot is tumour size. History suggests that
57
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.4: Kaplan-Meier Curves of Tumour Size (cm) categorised into 3 groups
0.0
Low to 1.5cm
1.51 to 2.5cm
2.51cm+
0
50
100
150
200
t
tumour size certainly is a factor in determining survival, so some time evaluating different
ways of using it to predict survival would be well spent. A variable of intrigue which
was touched on earlier was the tumour type, no emphasis had been put on this variable
as a potential prognostic factor but if an object has enough variation biologically to be
placed in to categories, then will its behaviour against survival still have enough variation
to be categorised? Eight of these types all cross over and follow very similar patterns of
survival, however there are two types which stand out. One has visibly better survival
time then the other types, but there are so few cases of this that we cannot statistically
prove any significant difference. And secondly there is a type which has visibly lower
survival times than the others. Now, this variable has more points, so there may be a
chance that this type is significantly different to the others. The method which would
usually be used to deduce this plays a large role in Survival Analysis and is known as the
Log-Rank Test.
Interesting finding: While age does not seem to play an important part in determining
survival based on Proportional Hazards (p=0.21); it may be of interest to see how age
58
0.0
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.5: Kaplan-Meier Curves of different types of Tumour
0
50
100
150
t
59
200
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.6: Kaplan-Meier Curves for different age groups
0.0
Under 40
In their 40's
In their 50's
60 and over
0
50
100
150
200
t
behaves as an ordinal categorical variable as opposed to a continuous one. Splitting
the age at first operation up into 4 groups: ≤ 40 (185 patients), 40 < to ≤ 50 (554
patients), 50 < to≤ 60 (761 patients) and > 60 (703 patients). As shall be seen later,
no significance was found in R through Cox Proportional Hazards(p=0.34). Every so
often however, intrigue persists, and the exploration of this idea continued. Figure 3.6
shows ploted survival estimates of the 4 age categories formed. Immediately, there was
graphical evidence that those diagnosed with breast cancer at younger than 40 would
have significantly less percentage survival than other age groups (especially between 5-8
years where the survival is a full 10% less than the other groups, which themselves do
not differ between each other in any significant way). A possible reason for this is that
more people at a younger age left the study (censored = 2) or another reason could be
the smaller numbers in that group. Kaplan-Meier estimates are invaluable in modern day
survival analysis, but there are a few restrictions which can occasionally make it difficult
to use them effectively:
1. Requires categorical predictors, as we saw earlier we had to order the continuous
variables into groups which can often lead to significant loss of information
60
2. Cannot accommodate time dependent variables; the cumulative survival approach
the method takes means that variables which vary over time will not be explained.
Thankfully the variables we are looking at are independent of time, so this is not
an issue
3. The don’t control for covariates
4. Descriptive, they are very useful at providing insight on the data but no real statistically significance can be formally derived
3.6.3
Life Table Method
The life table method (also known as the actuarial or Cutler Ederer method) is an approximation of the Kaplan-Meier method. It is based on grouped survival times and is
suitable for large data sets. The life table method assumes that subjects are withdrawn
randomly throughout each interval therefore, on average they are withdrawn half way
through the interval. This is not an important issue when the time intervals are short,
but bias may introduced when time intervals are long. This method also assumes that
the rate of failure within an interval is the same for all subjects and is independent of
the probability of survival at other time periods. Life tables are produced from large
scale population surveys (e.g. death registers) and are used less-frequently these days.
Generally speaking, the Kaplan-Meier method is preferred because it is less prone to bias.
3.6.4
Nelson-Aalen Estimators
The Nelson - Aalen estimator can be used to provide a non-parametric estimate of the
cumulative hazard rate function. Instantaneous hazard is dened as the proportion of the
population present at time t that fail per unit time. The cumulative hazard at time t, H(t)
is the summed hazard for all time up to time t. The relationship between cumulative
hazard and survival is as follows: H(t) = ln[S(t)] or S(t) = eH(t) The Nelson-Aalen
estimator of cumulative hazard at time t is dened as:
X
H(t) =
dj /rj , for 0 ≤ t ≤ t+
j:tj ≤t
Intuitively, this expression is estimating the hazard at each distinct time of death tj as the
ratio of the number of deaths to the number exposed.The cumulative hazard up to time t
is simply the sum of the hazards at all death times up to t, and has a nice interpretation
as the expected number of deaths in (0, t] per unit at risk. This estimator has a strong
justification in terms of the theory of counting processes.
61
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.7: Kaplan-Meier Curves for Endocrine Therapy
0.0
With Therapy
Without Therapy
0
50
100
150
200
t
3.7
Log-Rank Test
Previously, this report was exploring ways of clarifying whether there were significant
differences between groups in terms of survival. There was often strong suggestion from
the Kaplan-Meier curves that there were significant difference, but there was no way of
proving this using a simple plot. The Log-Rank test can be used perform comparisons
of survival curves using hypothesis tests. It is always of interest to measure the effect a
kind of therapy has on patients survival; The Kaplan-Meier plot in 3.7 shows that there
are obvious differences in the survival of patients who have had Endocrine (Hormonal)
therapy and those who have not. But the Log-Rank can be used to formally test if there
is significant between different groups, Rodriguez (2001). In strange way, making any
inferences of the performances of the therapies going on here is unlikely. These are not
Clinical trials, so are not comparing people who are not on the therapy against people in
the same condition who are not on the treatment (placebo). The patients who are not on
one kind of therapy will probably be on a range of other therapies. It is very difficult to
compare treatments for something as severe as cancer as there are many other underlying
factors such as the unpredictability of the tumour and other therapies the patient is
on. For example the situation in 3.7 had the people who were not on endocrine therapy
62
surviving longer than those who were on it, one might be thinking: Is this treatment really
working? In reality it is difficult to gauge a performance of this therapy as those who
are not using hormonal therapy may have further developed Cancer and using harsher
treatments such as Chemotherapy. Despite having a more intense treatment, their survival
chances will be less than those who do not require that treatment at that moment in time.
While its application in the above example is not in very useful for the purposes of our
study, there were some issues which were left earlier on which can be resolved using this
method. One of these was determining if there was significant difference in groups of
tumour type and also between the groups we created for the size of tumour.
Using the ”survdiff” function from the R survival package, the following output ensues
from Tumour Type:
Call:
survdiff(formula = surv ~ type, data = data)
n=2195, 8 observations deleted due to missingness.
N Observed Expected (O-E)^2/E (O-E)^2/V
type=1 1286
346
267.91
22.759
51.054
type=3
148
25
34.32
2.530
2.729
type=4
29
2
6.86
3.444
3.502
type=5
33
3
7.39
2.609
2.656
type=6
4
0
1.03
1.028
1.033
type=7
120
1
30.29
28.323
30.304
type=8
5
0
1.46
1.455
1.463
type=11 536
103
128.39
5.021
6.854
type=12
34
5
7.35
0.753
0.766
Chisq= 68.2
on 8 degrees of freedom, p= 1.12e-11
The p−value* shows straight away that there is significant difference between the groups.
Looking back at Figure 3.6, the younger age group seemed to have a slightly worse survival
time than the others - but at that time there was no concrete way to test these suspicions.
Step in the Log-Rank test (and R!)
* The p-value is the probability of obtaining a result as extreme or more extreme than the
one observed based on the current sample, given the null hypothesis is true.
Call:
survdiff(formula = surv ~ age.cat)
63
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure 3.8: Kaplan-Meier Curves for Menopausal Status
0.0
Pre−Menopause
Post−Menopause
0
50
100
150
200
t
age.cat=1
age.cat=2
age.cat=3
age.cat=4
N Observed Expected (O-E)^2/E (O-E)^2/V
185
54
39.6
5.243
5.723
554
123
128.1
0.206
0.281
761
155
168.4
1.063
1.632
703
153
148.9
0.113
0.164
Chisq= 6.6
on 3 degrees of freedom, p= 0.0842
As p= 0.0842 > 0.05, there is marginally insignificant difference in survival between these
age groups. That is not to say that there is NO difference, and there may well be other
factors contributing to these differences such as censoring. Menopause causes a great deal
of debate amongst scientists n breast cancer research. The original NPI did not consider it
an import factor, as have many other since but looking at the K-M plot 3.8, it seems like
there may be slight contrasts in survival between those who have/havent gone through
menopause.
Call:
survdiff(formula = surv ~ meno, data = data)
64
N Observed Expected (O-E)^2/E (O-E)^2/V
meno=1 817
190
185
0.1448
0.235
meno=2 1386
295
300
0.0892
0.235
Chisq= 0.2
on 1 degrees of freedom, p= 0.628
The log rank test above shows that there are no significant differences between the survival
of the two groups. Despite this, the menopausal status of a patient is still something
which could be looked into - it is obviously an indicator of age, but we have already seen
that age is not a large factor in survival. The K-M plot 3.8 shows that while there is
substantial cross over between the two groups, it is post-menopausal women who tend
to have the slight edge over those younger. Could this be something to do with the
fluctuating hormonal levels of a Pre-menopausal women? Perhaps there are some women
who have more volatile menstrual cycles, which can affect a women’s Oestrogen levels:
See the Appendix A.3. However, it should be remembered that Kaplan-Meier curves
are only informative and do not provide solid evidence of significance. Looking back
to the geography of Breast Cancer, it seemed that those who gave birth at an early
age in countries such as Mongolia were much less susceptible to breast cancer. These
are all interesting trends all of which could be relevant to breast cancer survival, and
hopefully continued research will reveal new information on Pre-Menopausal women and
their association with breast cancer.
65
Chapter 4
Cox Proportional Hazard Models
Cox Proportional Hazard Models are a class of survival models in statistics, first proposed
by Sir David Cox. Survival models relate the time that passes before some event occurs
to one or more covariates that may be associated with that quantity. In a proportional
hazards model, the unique effect of a unit increase in a covariate is multiplicative with
respect to the hazard rate. For example, taking a drug may halve one’s hazard rate for
a stroke occurring, or, changing the material from which a manufactured component is
constructed may double its hazard rate for failure. Other types of survival models such as
accelerated failure time models do not exhibit proportional hazards. These models could
describe a situation such as a drug that reduces a subject’s immediate risk of having a
stroke, but where there is no reduction in the hazard rate after one year for subjects who
do not have a stroke in the first year of analysis.
As mentioned, survival analysis typically examines the relationship of the survival distribution to covariates. Most commonly, this examination entails a specification of a
linear-like model for the log-hazard. For example, a parametric model based on the exponential distribution may be written as
loghi (t) = α + β1 xi1 + β2 xi2 + ..... + βk xik
(4.1)
or, equivalently
hi (t) = exp(α + β1 xi1 + β2 xi2 + ..... + βk xik )
(4.2)
that is, as a linear model for the log-hazard or as a multiplicative model for the hazard.
Here, i is a subscript for the observations, while the x’s are the co-variates. The constant
n this model represents a kind of log-baseline hazard, since loghi (t) = α when all of the
x’s are equal to zero.
66
The Cox model, in contrast leaves the baseline hazard function α(t) = logh0 (t) unspecified:
loghi (t) = α(t) + β1 xi1 + β2 xi2 + ..... + βk xik
(4.3)
or, again equivalently,
hi (t) = h0 (t)exp(β1 xi1 + β2 xi2 + ..... + βk xik )
(4.4)
This model is semi parametric because while the baseline hazard rate can take any form,
the covariates enter the model linearly. Now, if looking at two different observations (or
patients), say, i and i0 that differ in their x-values, with corresponding linear predictors
ηi = β1 xi1 + β2 xi2 + ..... + βk xik
(4.5)
and also
ηi0 = β1 xi0 1 + β2 xi0 2 + ..... + βk xi0 k
(4.6)
The Hazard ratio for these two observations,
hi (t)
h0 (t)eηi
=
hi0 (t)
h0 (t)eηi0
eηi
= η0
ei
ηi −ηi0
=e
(4.7)
which is independent of time t. Consequently, the Cox model is a proportional hazards
model. Even though the baseline hazard is unspecified, the cox model can still be estimated by the method of partial likelihood, developed by Cox in the same paper in which
he developed his original Proportional Hazards model.
4.0.1
Log Likelihood test
The Log Likelihood test will play a major part in the model building process. Its main
function is to compare two models, using hypothesis testing. The main criterion for this
test is that one model is nested in the other, with the nested model being the null hypothesis. An example would be comparing a model with 3 variables a, b, c against its
67
nested model with just two variables a, b to determine whether or not a models performance improves with the inclusion of the variable c. A test would then be performed to
determine whether or not the model is significantly better than the nested model. The
Log-Likelihood for the nested model is always going to be smaller than the alternative
model, but the p-value resting from the test will determine whether this significant or
not: If p < 0.05 then the null hypothesis is rejected and the alternative model is deemed
to be significantly better than the null model, conversely if p > 0.05 then the alternative
model is rejected so stick with the original model.
Very conveniently, this test is performed automatically when one is fitting models using
R; but to have s better idea of how it works, next is an example of the Log-Likelihood
test being used on the data.
Given a Proportional Hazards model with m variables
h(t | x1 , ...., xm ) = h0 (t)eβ1 x1 +β2 x2 +.....+βm xm
(4.8)
where xi are the co-variates such as Tumour Size, Menopausal Status etc. with βi their
respective coefficients and baseline hazard function h0 (t). The goal is to test whether or
not the model h(t | x) is improved with the inclusion of x1 =Tumour Grade.
Beginning with the null hypothesis
H0 : β1 = 0
(4.9)
H1 : β1 6= 0
(4.10)
With composite Hypothesis
This is essentially testing whether a proposed model
h(t | x1 ) = h0 (t)eβ1 x1
(4.11)
is significantly (α = 0.05) at least as accurate at predicting survival than the null model.
Using the test statistic
L(x0 | t, β)
D
= −2ln
L(x1 | t, β)
= −2ln(Likelihood of Null Model) + 2ln(Likelihood of Composite Model)
= 182.2
(4.12)
68
Reject H0 if D > χ21 at α = 0.05. Since comparing the difference of one parameter in the
model, there is 1 degree of freedom for the χ2 statistic, which is equal to 3.84 with 95%
confidence.
So D > 3.84 and reject H0 , concluding that x1 =Tumour Grade is an important factor in
determining survival.
4.0.2
Other Methods of Model Checking
Two alternative approaches for testing models is the Wald and Score statistics (which are
once again provided in the output of a model summary in R), which both also use the
likelihood to assess the models fit
Wald Statistic.... Again, say testing for β1 6= 0. The Wald test is that if
{
β1 2
} > χ21
SEβ1
(4.13)
then include x1 in the model.
This is far simpler than the above test, but for data that produce large estimates of the
coefficient, the standard error is often inflated, resulting in a lower Wald statistic, and
therefore the explanatory variable may be incorrectly assumed to be unimportant in the
model.
Score Statistic.... uses a simple null hypothesis that β is equal to some particular value
β1 .The statistic to test H0 : β = β1 is
S(β1 ) =
U (β1 )
I(β1 )
(4.14)
(4.15)
where
U (β1 ) =
∂logL(β1 | x)
∂β1
(4.16)
and fishers information
I(β1 ) = −E
∂2
logf (X; β1 ) | β1
∂ 2 β1
(4.17)
The score test once again uses S(β1 ) v χ21 to test the hypothesis. The main advantage
of the Score-test is that it does not require an estimate of the information under the
alternative hypothesis or unconstrained maximum likelihood.
69
4.1
Model Building using R
With this recently established understanding of the Proportional Hazards Model and how
it works, its now time to impute covariates in to the model in R and determine their
effectiveness by comparing them against the null model using the log likelihood ratio test.
The results are compared against a chi squared distribution at 95% confidence. This is
3.84 for 1 degree of freedom.
NB Care needs to be taken to not to confuse the log likelihoods when comparing models
in R; the output from a model automatically shows the log likelihood ratio against the
GLOBAL NULL model. The output will provide p-values for the test of each co-variate
being included/exclude in the model.
Three paths......remembering back, Histological Grade is made up of three components;
Mitosis, Tubular Formation and Pleomorphism. Perhaps one of these sub-variables is
more influential in predicting survival and can incorporated into the model. Figure 4.1
provides some insight. Tubular formation and Mitosis have very few patients of grade 1.
Grades 2 and 3 have good separation in Mitosis and Pleomorphism, less so with Tubular
formation, although the log-rank test shows significant difference (p= 4.29e-13). It might
be of interest to have a look at the survival of the summation of these grades in Figure
4.2, which shows the 7 groups. Remember the Tumour grade used in the current index
is calculated by the summation of these three factors, but can only take three values, so
perhaps a new variable of this form may maintain more information on the tumour. So
initially the model building process will be ran three times:
1. with the Histological Grade as a complete variable
2. with the Histological Grade removed and replaced by its components as three new
variables (M,P,T)
3. with the summation of the three component values (P+M+T)as a single variable
4.1.1
Model Building
The following shows the results of the full model (3) being run through the data in R.
The Variables associated with these R names can be seen in Table 1.8.
Call:
coxph(formula = surv ~ lymphnode + size + type + ER + ERc + radio +
axradio + chemo + endo + endotype + Op.type + detection +
axillary + total + LVI + age + nodespositive, data = data)
n= 2202, number of events= 485
70
Figure 4.1: Kaplan-Meier Curves of Histological Grades of Tumour, as well as of the
three factors involved in calculating the grade: Pleomorphism P, MitosisM and Tubule
Formation T
100
0.0
150
200
0
100
150
200
t
t
T
Histological Grade
0.0
0.4
^
S(t)
0.4
0.0
^
S(t)
50
0.8
50
0.8
0
0.4
^
S(t)
0.4
Grade I
Grade II
Grade III
0.0
^
S(t)
0.8
M
0.8
P
0
50
100
150
200
0
t
50
100
150
t
71
200
0.4
^
S(t)
0.6
0.8
1.0
Figure 4.2: Estimated Survival curves for 7 groups of new variable Total=T+P+M
0.0
0.2
T+M+P=3
=4
=5
=6
=7
=8
=9
0
50
100
150
t
72
200
coef exp(coef) se(coef)
z
lymphnode
0.321584 1.379311 0.102964 3.123
size
0.186466 1.204983 0.057591 3.238
type
0.019496 1.019687 0.013005 1.499
ER
-0.001122 0.998879 0.001072 -1.047
ERc
-0.105815 0.899591 0.185732 -0.570
radio
-0.104542 0.900737 0.077054 -1.357
axradio
0.081222 1.084611 0.068158 1.192
chemo
0.222189 1.248808 0.060984 3.643
endo
0.038994 1.039764 0.135264 0.288
endotype
-0.008773 0.991266 0.034322 -0.256
Op.type
0.112873 1.119490 0.107827 1.047
detection
0.280256 1.323469 0.121690 2.303
axillary
0.097829 1.102774 0.120183 0.814
total
0.346648 1.414319 0.044155 7.851
LVI
-0.413924 0.661051 0.103112 -4.014
age
0.012384 1.012461 0.005477 2.261
nodespositive 0.129774 1.138571 0.025323 5.125
--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1
lymphnode
size
type
ER
ERc
radio
axradio
chemo
endo
endotype
Op.type
detection
axillary
total
LVI
age
Pr(>|z|)
0.001788
0.001205
0.133856
0.295277
0.568869
0.174867
0.233389
0.000269
0.773131
0.798264
0.295190
0.021277
0.415646
4.11e-15
5.96e-05
0.023756
2.98e-07
1
exp(coef) exp(-coef) lower .95 upper .95
1.3793
0.7250
1.1272
1.6877
1.2050
0.8299
1.0764
1.3490
1.0197
0.9807
0.9940
1.0460
0.9989
1.0011
0.9968
1.0010
0.8996
1.1116
0.6251
1.2946
0.9007
1.1102
0.7745
1.0476
1.0846
0.9220
0.9490
1.2396
1.2488
0.8008
1.1081
1.4074
1.0398
0.9618
0.7976
1.3554
0.9913
1.0088
0.9268
1.0602
1.1195
0.8933
0.9062
1.3829
1.3235
0.7556
1.0426
1.6800
1.1028
0.9068
0.8713
1.3957
1.4143
0.7071
1.2971
1.5422
0.6611
1.5127
0.5401
0.8091
1.0125
0.9877
1.0017
1.0234
73
**
**
***
*
***
***
*
***
nodespositive
1.1386
0.8783
1.0834
Concordance= 0.77 (se = 0.014 )
Rsquare= 0.182
(max possible= 0.961 )
Likelihood ratio test= 442.5 on 17 df,
Wald test
= 478.5 on 17 df,
Score (logrank) test = 614.4 on 17 df,
1.1965
p=0
p=0
p=0
After fitting the three full models, it is evident that all are very similar so there is not
much to go on but initial model construction shall begin with choice (3) with total
(T+P+M). The Appendix contains the summary of the R output of the other two
models. The likelihood ratio test shows that models (1), (2) and (3) have likelihood ratio
scores of 433.6, 445.1 and 442.5 respectively. This is very little to go on at this point,
and on first glance the second model would be chosen. But remembering that this has
two extra variables, and operating on 19 df, so χ219 = 30.14 > χ217 = 27.59 which pushes
the choice into the favour of (3). In addition the second model has Pleomorpism(P) as
insignificant, to remove it at this early stage might be slightly ruthless.
LVI = x7 , chemo = x4 and total = x6 are all highly significant(***). lymphnode= x1 and
size = x2 are significant(**) terms while age = x8 and detection = x5 are marginally
significant(*). While Tumour Type = x3 is not significant at this point, it is not far off
and the analysis on it earlier suggest that it might yet play a role. So the first model with
8 co-variates is:
h(t | xi ) = h0 (t)exp(β1 x1 + β2 x2 + β3 x3 + β4 x4 + β5 x5 + β6 x6 + β7 x7 + β8 x8 )
(4.18)
Based on earlier analysis, Tumour Type 1 had significantly lower survival than the other
types, so type shall be re-categorise to be new.type where (1=Tumour Type 1, 2= Any
other tumour type). The same will be done for age and will be recategorised into a new
variable age.cat with 4 distinct categories 40 >,40 − 50,> 50 − 60 and 60 <.
Call:
coxph(formula = surv ~ lymphnode + size + new.type + chemo +
detection + total + LVI + age.cat, data = data)
n= 2195, number of events= 485
lymphnode
coef exp(coef) se(coef)
0.65777
1.93048 0.06761
74
z Pr(>|z|)
9.729 < 2e-16 ***
size
0.19650
1.21714 0.05515 3.563 0.000367
new.type
0.01430
1.01440 0.11937 0.120 0.904648
chemo
0.18015
1.19740 0.02856 6.307 2.84e-10
detection 0.33580
1.39906 0.12175 2.758 0.005812
total
0.35735
1.42954 0.04110 8.696 < 2e-16
LVI
-0.34246
0.71002 0.09974 -3.434 0.000596
age.cat
0.14389
1.15476 0.05005 2.875 0.004043
--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1
1
***
***
**
***
***
**
exp(coef) exp(-coef) lower .95 upper .95
lymphnode
1.930
0.5180
1.6909
2.2040
size
1.217
0.8216
1.0924
1.3561
new.type
1.014
0.9858
0.8028
1.2818
chemo
1.197
0.8351
1.1322
1.2663
detection
1.399
0.7148
1.1021
1.7761
total
1.430
0.6995
1.3189
1.5494
LVI
0.710
1.4084
0.5839
0.8633
age.cat
1.155
0.8660
1.0469
1.2738
Concordance= 0.761 (se = 0.014 )
Rsquare= 0.168
(max possible= 0.961 )
Likelihood ratio test= 404.3 on 8 df,
Wald test
= 390.9 on 8 df,
Score (logrank) test = 441.6 on 8 df,
p=0
p=0
p=0
lymphnode has now become the most significant variable; this is because there were other
variables which followed similar patterns and their removal has increased Lymph Nodes
power in the model. The re-categorised variable age.cat has seen it become significant
rather than marginally significant while new.type is still insignificant so it shall be removed from the model. LVI is the only negative term in the model, as it is numbered
1=Yes and 2=No; unlike all the other variables, if the value of LVI increases so does
the survival time for that patient. This is quite confusing, but it is merely a case of the
renumbering the variable. It was decided that chemo would be removed from the model.
Yes, it is a significant factor in determining survival, as one would expect: A patient who
undergoes Chemotherapy is likely to be in worse condition than someone who does nt,
thus it is a good indicator of survival. But this whole process is about establishing a
Prognostic Index which can weigh up all the current information on a patient soon after
they are diagnosed, and it is unlikely that they will have had any kind of Chemotherapy
by that point. Ideally a model would be constructed based around the information that
75
will be readily available for a patient as soon as possible so an Index score can be determined and the correct course of action can be taken. It does not make sense to be waiting
to see whether or no a patient is having Chemotherapy before telling them their survival
prospects.
The status of age.cat has also been reconsidered. Having been marginally insignificant
without categorisation, based on the exploratory analysis it was never going to be largely
significant in predicting survival. Yet, at this late stage, there it sits pretty in the model as
highly significant (***). This is not true, and after running the log-rank test on this new
variable it came to light that there were no significant differences between the 4 groups
(p=0.0842). So the second model is;
h(t | xi ) = h0 (t)exp(β1 x1 + β2 x2 + β5 x5 + β6 x6 + β7 x7 )
(4.19)
Call:
coxph(formula = surv ~ lymphnode + size + +detection + total +
LVI, data = data)
n= 2203, number of events= 485
coef exp(coef) se(coef)
z Pr(>|z|)
lymphnode 0.63653
1.88991 0.06813 9.343 < 2e-16
size
0.17292
1.18877 0.05481 3.155 0.001607
detection 0.18522
1.20349 0.11644 1.591 0.111660
total
0.33184
1.39353 0.03488 9.514 < 2e-16
LVI
-0.35262
0.70285 0.09856 -3.578 0.000347
--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1
1
exp(coef) exp(-coef) lower .95 upper .95
lymphnode
1.8899
0.5291
1.6537
2.1599
size
1.1888
0.8412
1.0677
1.3236
detection
1.2035
0.8309
0.9579
1.5120
total
1.3935
0.7176
1.3014
1.4921
LVI
0.7028
1.4228
0.5794
0.8526
Concordance= 0.754 (se = 0.014 )
Rsquare= 0.152
(max possible= 0.961 )
Likelihood ratio test= 362.7 on 5 df,
p=0
Wald test
= 341 on 5 df,
p=0
76
***
**
***
***
Score (logrank) test = 379.4
on 5 df,
p=0
Now, immediately with the removal of age.cat, the detection goes from being highly
significant to being not significant at all. This could perhaps be explained by a possible
relationship between older people going for screening more often than younger. In some
ways this makes sense, as regardless of whether a person has been screened or not, their
survival is always going to be determined by the progression and nature of the tumour.
The fact that someone has had their tumour detected earlier on through screening may
well reduce the potential risks which the factors can have, but any such changes will be
explained by the variables themselves. So there is no need to include method of detection
as all the information it can bring is already provided through the 4 remaining variables.
Call:
coxph(formula = surv ~ lymphnode + size + total + LVI, data = data)
n= 2203, number of events= 485
coef exp(coef) se(coef)
z Pr(>|z|)
lymphnode 0.63797
1.89264 0.06815 9.361 < 2e-16
size
0.18506
1.20330 0.05407 3.423 0.000620
total
0.34047
1.40561 0.03455 9.855 < 2e-16
LVI
-0.36264
0.69584 0.09856 -3.679 0.000234
--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1
1
***
***
***
***
exp(coef) exp(-coef) lower .95 upper .95
lymphnode
1.8926
0.5284
1.6560
2.1631
size
1.2033
0.8311
1.0823
1.3378
total
1.4056
0.7114
1.3136
1.5041
LVI
0.6958
1.4371
0.5736
0.8441
Concordance= 0.752 (se = 0.014 )
Rsquare= 0.151
(max possible= 0.961 )
Likelihood ratio test= 360.1 on 4 df,
Wald test
= 341.6 on 4 df,
Score (logrank) test = 378.9 on 4 df,
p=0
p=0
p=0
This leaves the final 4 factors. All of these are highly significant with the model, 3 of
which make up the original NPI. Notice that unlike the others, the coefficient for LVI is
negative. This is due to the fact that the variable was imputed with outcomes 1=Yes and
2=No, and moving in a positive direction numerically moves in a positive way in terms
77
of survival (No LVI −→ Better Survival), whereas as the other three variables move in a
positive direction, the survival moves in a negative direction. Which leaves the third and
final model:
h(t | xi )
= h0 (t)exp(β1 x1 + β2 x2 + β6 x6 + β7 x7 )
= h0 (t)exp(0.64x1 + 0.19x2 + 0.34x6 − 0.36x7 )
(4.20)
4.2
Index
Using the coefficients from the final model, an Index can be derived.
I = 0.64x1 + 0.19x2 + 0.34x6 − 0.36x7
(4.21)
Where
x1
= Lymph Node Involvement
x2
= Tumour Size (cm)
x6 = Mitosis+Pleomorphism+Tubular
x7
= Lymphovascular Invasion
(4.22)
4.3
Model Checking
After fitting an appropriate PH model for the given data, it is vital to check the goodnessof-fit and the residuals of the fitted model. The goodness-of-fit of a PH model mainly focuses on checking the validity of the assumption of the proportional hazards (i.e. whether
the effects of covariates on risk remain constant over time). Violations of PH assumption
can be resolved by:
• Stratification
• Adding time*covariate interaction
• Adding other time-dependent version of the covariate
4.3.1
Schoenfeld Residuals
For a more general Coxs regression model, the Proportional Hazards property is one of
the restrictions to using the model with time-fixed covariates. It assumes that the hazard
78
ratio between two sets of covariates is constant over time, because the common baseline
hazard function cancels out in the ratio of the two hazards. However, the impact of
time-varying covariates leads the hazard to vary over time, thus violating the assumption
of PH in the model. To test this PH assumption of Coxs regression model, Schoenfeld
(1980) introduced a dummy time-dependent covariate.
A Schoenfeld Residual......
is defined as the covariate value for the individual that failed minus its expected values
and yields residuals for each person that failed for each variable
OR
is the expected value of each covariate at time ti = a weighted average of the covariate,
weighted by the likelihood of failure for each individual in the risk set R at time tj
Schoenfeld residual
= xik −
Pj∈R(ti )
i=1
xjk pj
(4.23)
pj = Probability of event now for ith person
R(ti )
= Risk set
K
= Covariate of Interest
Xik
= Covariate-value for the person i
Plots of these residuals are useful in detecting non proportionality of predicted hazards
of the fitted model over the covariate space for each covariate. Under the proportional
hazards assumption, the Schoenfeld residuals have the sample path of a random walk;
therefore, they are useful in assessing time trend or lack of proportionality. Due to time
dependent covariates the generalized linear regression of the Schoenfeld residuals on functions of time gives a non-zero slope.Thus, a non-zero slope is an indication of a violation of
the proportional hazard assumption. As with any regression it is recommended to look at
the graph of the regression in addition to performing the tests of non-zero slopes. There
are certain types of non-proportionality that will not be detected by the tests of non-zero
slopes alone but that might become obvious when looking at the graphs of the residuals
such as nonlinear relationship (i.e., a quadratic fit) between the residuals.
In R, the cox.zph function (also found in the survival package) calculates tests of the
proportional hazards assumption for each covariate, by correlating the corresponding set
of Schoenfeld residuals with an appropriate transformation of time.
rho
chisq
p
lymphnode -0.057036 1.77e+00 1.84e-01
79
19
52
100
−2
0
2
4
● ●
●
●
●●
●●
●
●
●● ● ●●
●
●●●● ●●● ●
●
● ●
●●
●●
●
●
●
● ●●
● ●●
●●
●
●
●
●●
●
●
●●
● ●
●●
●●
●
●●● ● ●
●
● ●
●
●●
●
● ●●
●
●
●●●
●
● ●●
●●●● ●
●●●
●●● ●
●●
●● ●
●
● ● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
● ●
●●
● ●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
● ●●●
●● ●
●●
● ●
●●●
●
●●●●
●
●
●
●
●●●
● ●●●
●
●
●●
●
●
●
●
● ●●
●
●
●
●●
●
●●●●
●
●
●
●●
●●●
●
●●
●●
●
●●
●●
●
●●
●●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●● ●
●
●
●
●
●
●
● ●●●●
●● ●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●● ● ●
●●
●
●● ●●
●
●●
●
●●
●●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●● ● ● ● ● ●
● ●
●
●● ●
●●
●
●
●
●
●
●
● ● ●
●
●
●●
● ● ●●
●
●●
●●
●
●
●
●●
●
●
●●● ●●
●
●●●●
●
●●
●
● ●● ●●
●●●
●●●●
●
●
●●●●
●●
●● ● ●●●● ●●●
● ● ●●● ●●
●
●
●
●●●
●
●
●
●
●
● ●●
●● ●
Beta(t) for size
4
2
0
−2
Beta(t) for lymphnode
Figure 4.3: Plots of scaled Schoenfeld residuals against transformed time for each covariate
in the final CPH model
180
●
●
●
●●
●
● ● ●
●
●
● ●●●
●●●●
● ●●
●
●●
● ●● ● ●●
●●
●
●● ● ●
● ● ●●
●
● ●
●●
● ●● ● ●
●●
●
● ●●●●
●●
● ●●●
●
●
●
●●● ●
●● ●
●
●●
● ●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●● ● ● ●
● ●● ●●
●●
●●
●
●
●
●●●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
● ●●●
●
●●●
●
●●●●●
●●●● ●
●●
●
●●
●● ●● ●
●●
●
●
●
●
●
●●●●●●●
●
●●●●
●
●
●
●
●
●●
●●●
●
●●
●
●●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●● ●
●● ●
●●
● ●●
●
●
●●
●
●
●
●
●
●
●
●●
●
●●●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●
●
●●
●
●●
●
●
● ●● ●
●
●
●
●●
●
●
●●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●●
●●●
●● ●
●
●
●
●●
● ●●
●
●●●
●
●●
●●●
●
●
●
●
●
●
●●
●●●●● ●
●
●●●
●
●
● ●●●●●
●●
●●
●
●●
●
● ●
●
●
●●
●●
●
●
●●
● ● ●● ● ●
●●
●
● ●●
●
●
●●
●
●
●
19
52
52
100
0.005465
-0.344318
-0.000182
NA
4
2
0
−4
180
●
●
●● ● ● ● ● ● ● ● ● ●
●●●●
●
●●
●● ● ● ● ●
●
●● ●
●
●●●
●● ●
●
●
● ● ●●
●
●●
●●●●●
●
●
● ●
●
●
● ● ● ●
● ●
●
●
●● ●●
●
●
●
●
●
●
●
●
●
●
●
● ● ●●
●
●● ●
●
●
●
● ●●
●●●●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●●
●
●●
●
●
●●● ●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●● ● ●●● ●●●
●●●●
●●●●
● ● ●●● ●
●
●●
●●●●● ● ●
● ●●●● ●●
●●
●
● ●●●
●
●
●
●
●
●
●●
●●●
●
● ●●
●
●●
●
●
●
●
●
●
● ●●
●
●●
●
●
●●● ●
●
●
●●●
●●
●●●
●
● ●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●●
●
●
●●
●●
●
●
●●
●
●
●
●
●
●●●
●
●
●●●
●●
●
●●
●
●
●
●●
●
●
●
●
●
●●●
●●
●●
●●
●
●
●
●●
●
●
●
●
●
● ●●
●●●
●
●
●
●
●
●● ●
●●●●●
●●
●
●
● ●●
●● ●
●● ●
●
●
●●
●
●●● ●
●
●●●●
●●●●●●●
●●●
●●
●●
●●
●
●
●●●
●
●
●
●
●
●
●
●
●●●● ●●●
●
●
● ●
●
●● ● ●
●
19
Time
size
total
LVI
GLOBAL
−2
●
●●●
● ●
●●
● ●●●●
●
●
●
●
●
●
●●
●●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●●●●
●
●
●
●
●
●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●●
●
●
● ●
●● ● ●
●
●
●●
●●
●●
●●
●●
●●●
●
●●● ●
●
●●
●
●●●
●
●
●●
●
●
●
●●
●
●
●
●●
●
●●● ●
●
●●●●●
●
●
●●●
●●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●●● ●
● ●
●●●● ● ● ●
●
● ●
●
●
● ●●●● ● ●
● ●●
●
●●●●●●
●●
●
●
●
●● ●
● ●●● ●
●●
●●
●●
● ●●
●●● ● ●
●●●
●
●
●
● ●●
●
●
● ●
●
●
●●●●●
●
●
●● ●
●●●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●● ●● ●
● ● ●●●
●●
●●
● ●●
● ● ●●
●●
●● ●
● ●●
● ●●●●
●●
●
●●
●
●● ●●●●
●●
●●
●
●
● ●● ●
● ● ●
●● ●
●●
●●
●● ● ●
●
●
● ●
●
●
●
19
180
Time
Beta(t) for LVI
−0.5 0.5
−2.0
Beta(t) for total
1.5
Time
100
52
100
180
Time
1.41e-02
5.17e+01
1.67e-05
5.84e+01
9.06e-01
6.53e-13
9.97e-01
6.32e-12
The output R code above shows us that there is strong evidence of Non-proportional
hazards for the variable total. The object returned by cox.zph produces graphs of the
scaled Schoenfeld residuals against transformed time. Figure 4.3 shows the Schoenfeld
residuals, which are very useful in checking for violation of proportional hazards. The
solid line is a smoothing-spline fit to the plot, with the broken lines representing a ±2
standard error band around the fitted line. As expected from the original R output there
seems to be more of a trend in the residual plot of total than the other variables. One
method, Fox (2002) of perhaps attempting to remove this violation of proportional hazards
is to introduce an interaction term into the model. In this case, the model included an
interaction term between total and survival to try and remove the violation of PH.
As can be seen in the plot below this did not work out as planned, and now LVI and
lymphnode violate PH.
rho
chisq
p
80
lymphnode
0.09865
5.4777
size
0.02947
0.3935
total
0.00914
0.0807
LVI
-0.01023
0.0595
total:survival 0.31385 98.0970
GLOBAL
NA 455.8821
0.0193
0.5305
0.7763
0.8073
0.0000
0.0000
There is a lot of work going on by statisticians around the world in trying to refine
the model checking techniques for proportional hazards. The nature of these means
that more criteria have to be met than Linear and Generalized Linear models, such as
proportional hazards, which makes it difficult to construct a fully working model. For
a first attempt, having only one variable violating Proportional hazards is not bad. As
such, the prognostic index might well stand a chance of being a useful tool in predicting
survival. The interaction term will be swiftly removed and the final model is shown in
4.24. The higher the Index, the worse the Prognosis.
I = 0.64LN + 0.19T S + 0.34M P T − 0.36LV I
(4.24)
Where
LN
= Lymph Node Involvement
TS
= Tumour Size (cm)
M P T = Mitosis+Pleomorphism+Tubular
LV I
= Lymphovascular Invasion
81
Chapter 5
Summary, Conclusion & Further
Research
5.1
Summary and Conclusion
As with the all previous studies, Lymph Nodes was the most significant factor in determining survival (β1 = 0.64). The Histological Grade of the tumour was once again
a significant factor, although rather than be used represented as a three point categorisation of Tubular Formation (T), Mitosis (M) and Pleomorphism (P), it would instead
be represented as the 7-point categorisation of the summation (total in R), to prevent
loss of information. The 7-point categorisation seemed to work slightly better within the
model, with (β6 = 0.34 in the final Index. A third model was used which included T,P
& M each as individual variables, but it did not fair as well as the other two models.
One interesting finding did come from this model in that Pleomorphism was insignificant
within the model (p > 0.05) where as Mitosis and Tubular formation we very significant.
Lymphovascular Invasion (LVI) was the the fourth variable to come out significant after
the model building process, and was present in the final index with β7 = −0.36. The fact
that this variable is negative is purely down to the fact it was labelled different to the
other factors. Care should be taken if this Index is used on other datasets to make sure
that the co-variates are labelled the correct way around.
The recategorising of age into the variable(age.total) improved its predictive power
within the model, to a level of significance. Although this variable was eventually excluded
thorough lack of certainty over it, perhaps in the future more time can be exploring methods to include it in a prognostic index. The same can be said of a patients menopausal
status, the Kaplan-Meier estimates showed clear differentiation between pre-menopausal
and post-menopausal patients. The log-rank test rendered this difference marginally insignificant, but there is surely scope to find a way of incorporating these difference, in
however small a way, in a prognostic index. The fact that menopausal status is directly
82
correlated with age perhaps suggests the possibility that only one of them has a marginal
effect on survival, which is reflected in the other.
Further interesting points to come from exploratory analysis:
- The accuracy of measurement for tumour size seems questionable, with large numbers
represented for the factors of 0.5cm, specially for tumour sizes > 1.5cm. With tumour
size being a significant factor in the Nottingham Prognostic Index and in the Index
generated from this report (β2 = 0.19), this apparent rounding of tumour size effect the
5-year survival scores. It should be noted the knock on effect of this rounding will be
minimized in the NPI with Tumour Size weighted with a factor of 0.2 compared with 1
for Lymph node involvement and Histological grade
- There is a significant difference in survival for patients having their cancer detected
through screening against those who have their tumour detected through other methods.
Patients who have their cancer detected through screening generally identify the tumour
far earlier, hence improved survival. Figure 2.15 shows clear difference between the
survival % ’s, and Log-rank test shows clear difference between the groups (p < 0.05)
After deliberation, it was decided Method of detection should not be used in a prognostic
index because any decrease in the patients survival chances which result from late
detection would already be represented by the biological attributes of the tumour. Those
patients who do have regular screening have on average 10% better survival, yet only
32.8% have regular checks. This is perhaps the most striking statistic to come from this
report, and findings like this should act towards encouraging the female population to
attend regular check ups
As well as covering all of the main areas of survival analysis, this report has touched
upon less well known techniques. Kaplan-Meier curves proved the most useful method
of gaining insight into the behaviour by providing very readable visualisations of the
relationship of variables with survival. These visual plots are also beneficial in that
they can be used to explain findings to those less familiar with the statistical techniques
behind this research. Cox-Proportional Hazard models provide a very efficient way of
establishing significant variables, and using R one can very quickly single out important
pieces of information on their behaviour within the model. The correct approach to take
for checking a model remains an issue, although this report demonstrated their usefulness
in determining violations of Non proportionality.
5.2
Further Research
This has been a fascinating study to be a part of, and unfortunately there was simply
not enough time to achieve all that was hoped from the many points covered in this
83
project, each of which can open the door for further research. Such is the nature of these
modern statistical techniques, new questions are always arising and time would be well
spent seeking the answer to them.
1. The Prognostic Index designed (4.24) in Chapter 4 currently has no function other
than providing information on the factors involved with survival, and their comparative significance through their β-values. With further research, this Index could be
configured against data to separate the Index into groups and associate these with
% survival scores for different gradings constructed. With this Index and associated Survival Scores, this model could be validated by running patients though the
score and comparing their predicted survival with actual survival (and indeed that
obtained from the NPI). One point which was touched upon earlier in the report
was the possibility of rather than using 5-year scores, using a 7-year survival score.
This could also be explored
2. The four factors which were established as being important in predicting survival
came from a proportional Hazards model which violated the assumption of proportional hazards. More time could be invested in finding a way to augment the model
and the variables inside of it to improve its accuracy and goodness of fit. This report
explored Schoenfeld residuals as a method for testing a fit of a CPH model, (which
in themselves are many different areas to study) but there are different methods
out there such as Martingale residuals which could also be employed in the model
checking stages.
3. More of the same research looking at Centres other than (1), as was seen figure 1.8
there are quite significant differences in the average Nottingham Prognostic Index
scores. This could be for a number of reasons such as different methods of treatment,
or perhaps specialized for specific kinds of breast cancer, all of which would be very
interesting to look at.
4. Continue with the scrutiny of the NPI, and continue to gain insight into its shortfalls
and what can be learnt from it
5. Read further into the most effective way the three components of Histological Grade
can be used as predictive elements. Pleomorphism was insignificant when isolated as
a single variable, but more research could go into whether it is wholly insignificant
or if it is significant when used with other factors
6. Expand upon the interesting findings on the measurement of tumour size, and deduce whether this could truely be an issue in misrepresentation of a patient within
a prognostic index
84
5.3
Acknowledgements
I would like to thank Professor Andy Wood of Nottingham University for his continued
help and support throughout this eventful project. Special thanks must also go to Dr.
Steven Chan of Nottingham City Hospital and Professor Graham Balls for giving me the
opportunity to work on a dissertation of this nature.
85
Appendices
86
Appendix A
Some further Interesting findings
R Output of Full model (1) with Tumour Grade chosen.
Call:
coxph(formula = surv ~ lymphnode + size + grade + type + ER +
ERc + radio + axradio + chemo + endo + endotype + Op.type +
detection + axillary + LVI + age + nodespositive, data = data)
n= 2202, number of events= 485
(1 observation deleted due to missingness)
coef exp(coef)
lymphnode
0.299822 1.349618
size
0.197924 1.218870
grade
0.714215 2.042582
type
0.015806 1.015932
ER
-0.001346 0.998655
ERc
-0.155817 0.855716
radio
-0.105538 0.899840
axradio
0.083471 1.087053
chemo
0.221175 1.247542
endo
0.037318 1.038023
endotype
-0.012185 0.987889
Op.type
0.122273 1.130062
detection
0.291767 1.338791
axillary
0.092484 1.096896
LVI
-0.407942 0.665017
age
0.011634 1.011702
nodespositive 0.132179 1.141312
---
se(coef)
0.103543
0.056937
0.097810
0.012853
0.001070
0.184748
0.074691
0.065329
0.060331
0.135088
0.034601
0.107681
0.121626
0.120481
0.102926
0.005479
0.025645
87
z
2.896
3.476
7.302
1.230
-1.257
-0.843
-1.413
1.278
3.666
0.276
-0.352
1.136
2.399
0.768
-3.963
2.124
5.154
Pr(>|z|)
0.003784
0.000509
2.83e-13
0.218791
0.208654
0.399004
0.157659
0.201357
0.000246
0.782356
0.724723
0.256162
0.016445
0.442710
7.39e-05
0.033712
2.55e-07
**
***
***
***
*
***
*
***
Figure A.1: Histogram of Survival for patients dying within 10 years of prognosis with
predicted survival of 85%
10
120.00%
9
100.00%
8
Number of Patients
7
80.00%
6
5
60.00%
Frequency
4
Cumulative %
40.00%
3
2
20.00%
1
0
0.00%
NPI Scores
Signif. codes:
lymphnode
size
grade
type
ER
ERc
radio
axradio
chemo
endo
endotype
Op.type
detection
axillary
LVI
age
0 *** 0.001 ** 0.01 * 0.05 . 0.1
1
exp(coef) exp(-coef) lower .95 upper .95
1.3496
0.7410
1.1017
1.6533
1.2189
0.8204
1.0902
1.3628
2.0426
0.4896
1.6863
2.4742
1.0159
0.9843
0.9907
1.0419
0.9987
1.0013
0.9966
1.0008
0.8557
1.1686
0.5958
1.2291
0.8998
1.1113
0.7773
1.0417
1.0871
0.9199
0.9564
1.2355
1.2475
0.8016
1.1084
1.4041
1.0380
0.9634
0.7966
1.3527
0.9879
1.0123
0.9231
1.0572
1.1301
0.8849
0.9150
1.3956
1.3388
0.7469
1.0548
1.6992
1.0969
0.9117
0.8662
1.3891
0.6650
1.5037
0.5435
0.8137
1.0117
0.9884
1.0009
1.0226
88
0.0
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure A.2: Kaplan-Meier Curves of NPI with 95% CI
0
50
100
150
t
89
200
0.2
0.4
^
S(t)
0.6
0.8
1.0
Figure A.3: Kaplan-Meier estimates of Estrogen Receptor Status
0.0
ER Code Negative
ER Code Positive
0
50
100
150
t
90
200
nodespositive
1.1413
0.8762
1.0854
1.2001
Concordance= 0.772 (se = 0.014 )
Rsquare= 0.179
(max possible= 0.961 )
Likelihood ratio test= 433.6 on 17 df,
p=0
Wald test
= 467.1 on 17 df,
p=0
Score (logrank) test = 606 on 17 df,
p=0
R summary of full model (2) with individual variables T,P,M included
Call:
coxph(formula = surv ~ lymphnode + size + type + ER + ERc + radio +
axradio + chemo + endo + endotype + Op.type + detection +
axillary + T + P + M + LVI + age + nodespositive, data = data)
n= 2202, number of events= 485
(1 observation deleted due to missingness)
coef exp(coef) se(coef)
z
lymphnode
0.321164 1.378732 0.102950 3.120
size
0.181484 1.198996 0.057674 3.147
type
0.019126 1.019310 0.013274 1.441
ER
-0.001077 0.998924 0.001075 -1.002
ERc
-0.105176 0.900166 0.185757 -0.566
radio
-0.099991 0.904846 0.076269 -1.311
axradio
0.076435 1.079432 0.067495 1.132
chemo
0.221688 1.248182 0.061133 3.626
endo
0.045326 1.046369 0.135799 0.334
endotype
-0.006510 0.993511 0.034355 -0.189
Op.type
0.113824 1.120555 0.107772 1.056
detection
0.282151 1.325979 0.121720 2.318
axillary
0.088392 1.092417 0.120346 0.734
T
0.384513 1.468898 0.106985 3.594
P
0.154848 1.167480 0.126398 1.225
M
0.424319 1.528548 0.074003 5.734
LVI
-0.414804 0.660470 0.103185 -4.020
age
0.011928 1.012000 0.005483 2.176
nodespositive 0.130908 1.139863 0.025316 5.171
--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1
91
Pr(>|z|)
0.001811
0.001651
0.149618
0.316378
0.571255
0.189844
0.257448
0.000287
0.738552
0.849709
0.290897
0.020447
0.462651
0.000326
0.220545
9.82e-09
5.82e-05
0.029579
2.33e-07
1
**
**
***
*
***
***
***
*
***
exp(coef) exp(-coef) lower .95 upper .95
lymphnode
1.3787
0.7253
1.1268
1.6870
size
1.1990
0.8340
1.0708
1.3425
type
1.0193
0.9811
0.9931
1.0462
ER
0.9989
1.0011
0.9968
1.0010
ERc
0.9002
1.1109
0.6255
1.2955
radio
0.9048
1.1052
0.7792
1.0507
axradio
1.0794
0.9264
0.9457
1.2321
chemo
1.2482
0.8012
1.1072
1.4071
endo
1.0464
0.9557
0.8018
1.3655
endotype
0.9935
1.0065
0.9288
1.0627
Op.type
1.1206
0.8924
0.9072
1.3841
detection
1.3260
0.7542
1.0445
1.6832
axillary
1.0924
0.9154
0.8629
1.3830
T
1.4689
0.6808
1.1910
1.8116
P
1.1675
0.8565
0.9113
1.4957
M
1.5285
0.6542
1.3222
1.7671
LVI
0.6605
1.5141
0.5395
0.8085
age
1.0120
0.9881
1.0012
1.0229
nodespositive
1.1399
0.8773
1.0847
1.1978
Concordance= 0.771 (se = 0.014 )
Rsquare= 0.183
(max possible= 0.961 )
Likelihood ratio test= 445.1 on 19 df,
p=0
Wald test
= 483 on 19 df,
p=0
Score (logrank) test = 624.8 on 19 df,
p=0
92
Bibliography
[1] D.R.Cox & D.Oakes, 1984, Analysis of Survival Data, Chapman & Hall
[2] G. Rodrguez, 2001, Non-Parametric Estimation in Survival Models, Princeton Press
[3] J.Fox, 2002, Cox Proportional-Hazards Regression for Survival Data
[4] R.J.Little & D.B.Rubin, 2002, Statistical Analysis with Missing Data, WileyInterscience
[5] M.Galea,R.Blamey,C.Elston,I.O.Ellis, 1992,The Nottingham Prognostic Index in Primary breast Cancer Kluwer Academic publishers
[6] N.E Breslow, 1976, Analysis of Survival Data under the Proportional Hazards
Model,International Statistical Review
[7] Schoenfeld, D , 1981, The asymptotic properties of nonparametric tests for comparing
survival distributions, Biometrika
[8] ,Kaplan,E.L.& P.Meier, 1958. Non-parametric Estimation from Incomplete Observations. J. Am. Stat. Assoc., 53:457481
[9] Allison,P D, 2001 Missing Data (1st ed.), Thousand Oaks: Sage Publications,
[10] D. Mayer, 2004, Essential Evidence Based Medicine, Cambridge Press Publication
[11] R.B Millar, 2011, Maximum Likelihood Estimation and Inference: With Examples in
R, SAS and ADMB,John Wiley & Sons,
[12] Kendall, M.G., Stuart, A. 1973, The Advanced Theory of Statistics, Volume 2 Griffin.
[13] N. Breslow & J. Crowley, 1974, A Large Sample Study of the Life Table and Product
Limit Estimates Under Random CensorshipAnn. Statist. Volume 2, Number 3
93