Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2963459

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Cure rate modelling under latent activation schemes
- 3 Latent spatial cure rate model
- 4 Bayesian estimation of spatial cure models
- 5 Illustrations
- 6 Future work
- References

Authors

Related links

Stat Methods Med Res. Author manuscript; available in PMC 2010 October 25.

Published in final edited form as:

PMCID: PMC2963459

NIHMSID: NIHMS127620

Freda Cooner, Division of Biostatistics, School of Public Health, University of Minnesota, Minneapolis, MN, USA;

Address for correspondence: Sudipto Banerjee, Division of Biostatistics, School of Public Health, University of Minnesota, Mayo Mail Code 303, Minneapolis, Minnesota 55455-0392, USA. ude.nmu.tatsoib@botpidus

The publisher's final edited version of this article is available at Stat Methods Med Res

See other articles in PMC that cite the published article.

The emergence of geographical information systems and related softwares nowadays enables medical databases to incorporate the geographical information on patients, allowing studies in *spatial* associations. Public health administrators and researchers are often interested in detecting variation in survival patterns by region or county in order to understand the possible factors that contribute towards such spatial discrepancies. These issues have led statisticians to develop survival models that account for spatial clustering and variation. Additionally, with rapid developments in medical and health sciences, researchers increasingly encounter data sets where a substantial portion of patients are *cured*. Models accounting for cure in the population assist in the prognosis of potentially terminal diseases. This article proposes a Bayesian modelling framework that models spatial associations for areally referenced survival data using a general class of cure models proposed by Cooner *et al.* The special models we outline are alternatives to the traditional proportional hazards models and can be fitted using standard Bayesian software such as WinBUGS.

With the advent of geographical information systems, medical and health databases now often include geographical information on patients (e.g., county of residence) that allow the researchers to investigate *spatial* associations in health phenomena. Spatial disease mapping, an area that has attracted much attention over the last several years,^{1}^{–}^{3} concerns the development of statistical and cartographic methods that map disease rates accounting for different types of extraneous variation. Lawson and Williams,^{3} in particular, discuss several aspects of modelling and software implementations.

In recent times, there has been growing interest in capturing spatial trends in the *survival* patterns of patients suffering from potentially terminal diseases. Health care administrators and public health researchers are often interested in detecting variation in survival patterns by counties for determining the possible factors that contribute towards such variability. These issues have led statisticians to develop survival models that account for spatial clustering and variation. Banerjee *et al.*^{4} investigated the spatially correlated frailties in traditional parametric survival models adopting a hierarchical Bayesian approach. Instead, Li and Ryan^{5} address this problem from a classical perspective.

Generally, the spatial survival modelling mentioned above is treated in a proportional hazards framework. With rapid developments in medical and health sciences, scientists and health professionals encounter more data sets where the patients are expected to be *cured*. Formulation and estimation of models that account for cure are important for understanding prognosis in potentially terminal diseases. Traditional parametric survival models such as Weibull or Gamma^{6} do not account for cure, assuming instead that individuals who do not experience the event are *censored*. Although subtle, in these contexts one distinguishes between the concepts of censoring and cure: censoring refers to a subject who does not fail within the time window of the experiment, while cure refers to one who will *never* fail. Indeed the latter is an abstraction as we never ‘observe’ a cure (due to a finite monitoring time). Still estimating the probability of such an outcome, especially in various cancer-relapse settings, can help expose unknown health issues concerning that population.

Statistical models address this conceptually challenging problem by parametrizing the probability of a cure, called the *cure fraction*. One of the earliest such models was a class of mixture models by Berkson and Gage^{7} that generated subsequent investigations by Farewell,^{8}^{,}^{9} Goldman,^{10} and Ewell and Ibrahim^{11} among others. Yakovlev *et al.*^{12} and Chen *et al.*^{13} offered an alternative approach to formulating cure models^{14}^{,}^{15} that assume a latent biological process in generating the observed failure (say, cancer relapse). Recently, Cooner *et al.*^{16} proposed a very general class of cure models assuming (possibly several) *latent factors* or *latent risks* (e.g., metastasis-competent tumours that lead to cancer relapse) corresponding to each patient. For an individual to be at risk of failure, one must be exposed to at least one of these latent factors. If not, the individual is not at risk and is considered *cured*. Failure is observed when some (perhaps one or all) of these latent factors become *activated*.

In the modelling of spatial association using random effects and, more generally, regressors in cure models, the modeller faces some conceptually interesting alternatives. Certain data sets may be better modelled by assuming that the cure fractions are spatially associated. Alternatively, one may model regressors and spatial effects in the distribution of the latent factors. Richer models incorporating such random effects in both the cure fraction and the latent distributions can also be conceptualized. However, identifiability and estimability of such models raise concerns that have hitherto gone unaddressed with most of the spatial-survival analysis literature focussing upon proper survival models. Recently, Banerjee and Carlin^{17} demonstrated data analysis from a spatially referenced interval-censored smoking-cessation study in southern Minnesota using a binary cure model interpreted as *one* latent factor. However, their models preclude studying spatial effects in the cure fraction, as such effects are not estimable.

This article explores spatial cure rate models by extending the modelling framework of Cooner *et al.*^{16} We propose a methodological framework that enables versatile and flexible modelling of spatial associations for survival data in which a prominent proportion of subjects might be cured. We adopt a Bayesian hierarchical framework that allows very rich modelling structures estimated using Markov Chain Monte Carlo (MCMC) methods.^{18} It is important to point out that we limit our discussion to *areal* settings where the geographical referencing for each subject is by the county they live in, rather than by the coordinates of each subject’s residence (*point-referencing*). Areal summaries are frequently found in public health data due to data privacy issues.

With information at the county level, spatial modelling will proceed from the neighbourhood information of the counties. More precisely, we consider a vector of univariate latent random variables ϕ = (ϕ_{1},…,ϕ_{n}) in *n* areal locations (such as counties, states) that follow a *spatial distribution*. Unlike the independent random effects that specify ϕ ~ *N*(0, τ^{2}*I*), a spatial distribution will model ϕ ~ *N*(0, Σ^{−1}), where Σ^{−1} is a choice of the spatial precision matrix. Typically, these effects are modelled after adjusting for covariates or risk information and hence are zero-centered. Specifying Σ^{−1} based on the underlying neighbourhood configurations follow from the theory of Markov Random Fields (see e.g., the works of Rue and Held^{19} for details). We work with univariate and multivariate conditionally autoregressive (CAR) models that accrue computational benefits, while incorporating spatial effects. In particular, we investigate a Generalized multivariate CAR (GMCAR) distribution, recently investigated by Jin *et al.*^{20} that encompasses very rich multivariate spatial associations.

The remainder of this paper evolves as follows. In Sections 2 and 3 we outline the cure rate and spatial modelling framework that we work with. Section 4 discusses the Bayesian implementation algorithms. Section 5 provides the analysis of a breast cancer survival data set obtained from the National Cancer Institute’s Surveillance, Epidemiology and End Results (SEER) database (http://www.seer.cancer.gov/). Finally, Section 6 summarizes and indicates future areas for research.

Cure models based upon latent activation schemes^{16} involves failure times at two different levels: an *observed* failure time *T* corresponding to the *observed* time when an individual *fails*, as well as the *latent* event times, *Y _{k}, k* = 1,…,

To derive the distribution of *T*, we assume that *r* out of *N* latent factors need to be activated for the subject to fail, so *T* = *Y*_{(r)}, *r* = 1,…, *N* where *Y*_{(1)} < … < *Y*_{(N)} are the ordered *Y _{k}*’s. Cooner

The conditional distribution of *T* given *N* can be written down in terms of the incomplete beta function, or a beta cdf denoted by *IB*(*S*(*t*); *N* − *r* + 1, *r*), (using a standard result on order statistics using the binomial theorem, e.g., Rao^{21} (p. 215) as

$$P(T\ge t|N)=1(N=0)+IB(S(t);N-r+1,r)1(N\ge r\ge 1),$$

(1)

where 1(·) is the indicator function and

$$\begin{array}{ll}IB(S(t);N-r+1,r)\hfill & =\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{j=0}^{r-1}\left(\begin{array}{c}N\\ j\end{array}\right){[F(t)]}^{j}{[S(t)]}^{N-j}}\hfill \\ \hfill & =\phantom{\rule{thinmathspace}{0ex}}N\left(\begin{array}{c}N-1\\ r-1\end{array}\right){\displaystyle {\int}_{0}^{S(t)}{u}^{N-r}{(1-u)}^{r-1}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}u.}\hfill \end{array}$$

The unconditional survival function of *T*, *S**(*t*), is related to the latent distribution as

$${S}^{*}(t)={E}_{N}[P(T\ge t|N)]=P(N=0)+{E}_{N}[IB(S(t);N-r+1,r)1(N\ge r\ge 1)],$$

(2)

where the expectation is taken over the distribution of *N*. Here *S**(*t*) always exists, being bounded between 0 and 1 for any valid distribution of *N* restricted to *N* ≥ *r* ≥ 1 as *IB*(*S*(*t*); *N* − *r* + 1, *r*) is otherwise 0. Also, as lim_{t→∞} *S*(*t*) = 0, we have lim_{t→∞} *S**(*t*) = *P*(*N* = 0), showing that *S**(*t*) is *improper* whenever *P*(*N* = 0) > 0. Indeed, *P*(*N* = 0) is the probability of a person being cured or immune, hence called the *cure fraction*, and depends only upon the distribution of *N*, irrespective of what *r* is. Despite *S**(*t*) being improper, the hazard, say *h**(*t*)*dt* ≈ *P*(*T* [*t*, *t* + d*t*)|*T* > *t*), is still valid so that *h**(*t*) is evaluated as −(d/d*t*) log *S**(*t*), or *f**(*t*)/*S**(*t*), where *f**(*t*) is the corresponding *improper* density. In fact, if *f* (*t*) is the proper latent density corresponding to *F*(*t*), then

$${f}^{*}(t)=f(t){E}_{N}\left[N\left(\begin{array}{c}N-1\\ r-1\end{array}\right){[S(t)]}^{N-r}{[F(t)]}^{r-1}1(N\ge r\ge 1)\right].$$

(3)

The variable *N* can never be *observed* and must be modelled using a probabilistic assumption. From a more theoretical perspective, one can show^{16} that if *N* ~ Po(θ) and if *r* or *N* − *r* is fixed at some positive integer (e.g., *r* = 1 or *r* = *N*), then θ is identifiable under a scale-invariant prior *g*(θ) 1/θ permitting regression in log(θ). These generate a well-identified class of cure models, obtained by setting *r* appropriately. Note that with *N* ~ Po(θ), we have the cure fraction exp(−θ). This theoretical identifiability permits regression (along with spatial random effects) in log(θ), amounting to a Poisson regression.

The first-activation scheme (*r* = 1) assumes that activation of a single latent factor will lead to observed failure. According to a biological model for patients diagnosed with cancer, *N* is the number of metastasis-competent clonogenic cells that are in an irreversible process towards metastasis, and *Y _{k}* is the time for the

The physical framework of the first-activation scheme is not unique in assuming that one latent factor’s activation generates observed failure. Alternatively, *N* can be the number of latent factors that must *all* be activated for failure. For instance, biological models for certain types of cancer posit that a patient’s immune response gets activated after the initiation of *N* cell mutations. This immune response may be able to resist up to *N* − 1 promotions of mutated cells before disease manifestation or death. Here, we model *Y _{k}* as the time to promotion of the

Turning to alternative modelling of *N*, we note that several authors including Tucker *et al.*^{22} and Tucker and Taylor,^{23} have questioned the Poisson assumption for any particular cancer cure scenario. From a modelling standpoint, the number of possible latent events *N* can have any finite-mean integer-valued distribution (e.g., binary, geometric and so on). In traditional cure rate models,^{7}^{,}^{9} *N* is binary with only one latent dominant event (e.g., a dominant metastasis competent tissue-mass in breast cancer), so *N* ~ Ber(θ) (with θ being the probability of an activation). Here *S**(*t*) = 1 − θ(1 − *S*(*t*)), approaches the cure fraction 1 − θ.

Another specification we consider here is *N* ~ Geo(θ), a geometric distribution with mean θ/(1 − θ)(*P*(*N* = *n*) = θ^{n}(1 − θ)). A biological motivation, different from the clonogen motivation in ^{Chen et al. (1999)}^{,}^{13} can be offered as follows (also see the work of Moolgavkar *et al.*^{24}). Assuming a short time-interval of the mutation (initiation) period (due to exposure to genetic damage), the patient’s body produces a sequence of *N* mutated cells/tissues before activating the immune system. Every mutation (initiation) may give rise to an effective immune response from the body with probability 1 − θ, capable of destroying the last mutated tissue/cell and halting the mutation process. With ${\left\{{Y}_{k}\right\}}_{k=1}^{N}$ now being the promotion times of the mutated cells, a first-activation scheme with *T* = min_{k} *Y _{k}* models failure with any one activation. Now we have

The cure rate models in Cooner *et al.*^{16} accommodate two primary regression structures, one using the mean of the latent factors, viz. the *Y _{k}*’s, or using a Poisson regression in the cure fraction itself. More precisely, suppose we have

For studying spatial variation across regions and, more specifically, smoothing across regions using adjacency information, we introduce *spatial frailties*^{4} that are spatially correlated region-specific random effects ϕ_{i}, *i* = 1,…, *I*. Interesting modelling choices arise. We could, as for the regressors, add these frailties either in η or in log(θ). If geographical variation is expected to reflect itself in lurking covariates that affect the latent factors, we opt for the former. In contrast, if we expect this variation in the cure fraction itself, we opt for the latter. However, in practice this choice will often not be clear. Therefore, more generally, one may incorporate two different sets of frailties, ${\left\{{\varphi}_{1i}\right\}}_{i=1}^{I}$ in the η_{ij}’s and ${\left\{{\varphi}_{2i}\right\}}_{i=1}^{I}$ in the log(θ_{ij})’s. On the basis of these, we classify models as univariate or bivariate; we discuss them in detail below.

In principle, *N* can be modelled using any integer-valued probability distribution. Irrespective of the distribution of *N*, the regression coefficients and frailties are identifiable from the data when they model the η_{ij}, giving rise to the following models:

$$\begin{array}{c}\text{\hspace{1em}Model}\phantom{\rule{thinmathspace}{0ex}}1(\mathrm{a}):{\eta}_{ij}={\mathbf{x}}_{ij}^{\mathrm{T}}\beta +{\varphi}_{i};\phantom{\rule{thinmathspace}{0ex}}{\theta}_{ij}=\theta ;\text{}N~\text{Po}(\theta );\\ \text{\hspace{1em}Model}\phantom{\rule{thinmathspace}{0ex}}1(\mathrm{b}):{\eta}_{ij}={\mathbf{x}}_{ij}^{\mathrm{T}}\beta +{\varphi}_{i};\phantom{\rule{thinmathspace}{0ex}}{\theta}_{ij}=\theta ;\text{}N~\text{Ber}(\theta );\\ \text{\hspace{1em}Model}\phantom{\rule{thinmathspace}{0ex}}1(\mathrm{c}):{\eta}_{ij}={\mathbf{x}}_{ij}^{\mathrm{T}}\beta +{\varphi}_{i};\phantom{\rule{thinmathspace}{0ex}}{\theta}_{ij}=\theta ;\text{}N~\text{Geo}(\theta );\\ \text{Model}\phantom{\rule{thinmathspace}{0ex}}2:\phantom{\rule{thinmathspace}{0ex}}\text{log}({\theta}_{ij})={\mathbf{x}}_{ij}^{\mathrm{T}}\beta ;\phantom{\rule{thinmathspace}{0ex}}{\eta}_{ij}={\eta}_{0}+{\varphi}_{i};\text{}N~\text{Po}(\theta ).\end{array}$$

In contrast, covariates modelled through regression on the cure fraction will be identified only if *N* ~ Po(θ), thereby precluding models 1(b)–(c). The spatial cure rate models for studying the spatial associations in cure fractions are given by

$$\begin{array}{r}\hfill \text{Model}\phantom{\rule{thinmathspace}{0ex}}1(\mathrm{a}):{\eta}_{ij}={\mathbf{x}}_{ij}^{\mathrm{T}}\beta ;\phantom{\rule{thinmathspace}{0ex}}\text{log}({\theta}_{ij})={\theta}_{0}+{\varphi}_{i};\text{}N~\text{Po}(\theta );\\ \text{Model}\phantom{\rule{thinmathspace}{0ex}}2:\text{log}({\theta}_{ij})={\mathbf{x}}_{ij}^{\mathrm{T}}\beta +{\varphi}_{i};\phantom{\rule{thinmathspace}{0ex}}{\eta}_{ij}=\eta ;\text{}N~\text{Po}(\theta ).\hfill \end{array}$$

Note that the interpretation of the parameters will be different for the two sets of models. In models 1(a)–(c), higher values of η_{ij} indicate longer latent event times and hence prolonged survival, so covariates influence the latent factors that cause failure. In contrast, in model 2, covariates influence the cure fraction directly. Almost certainly, most covariates that affect the latent event times will also affect the cure. For these covariates, we have a consistency in behaviour under the two groups of models. To see this, note that the Weibull mean is a decreasing function of η_{ij} (given by exp(−η_{ij}/ρ)Γ(1 + 1/ρ)) and the cure fraction is a decreasing function of θ_{ij}. This is sensible because covariates that adversely affect the survival time will likely affect the cure fraction adversely as well (and vice versa). Therefore, it is expected that their coefficients have the same sign under both sets of models.

Turning to the spatial modelling of the ϕ_{i}’s, we employ variants of the (univariate) conditionally autoregressive model introduced by Besag.^{25} For further theoretical discussions, see the works of Cressie,^{26} Banerjee,^{27} Rue and Held,^{19} Heikkinen and Högmander,^{28} Högmander and Møller^{29} and Hoeting *et al.*^{30} for generalizations and applications. Collecting the frailties into an *I* × 1 vector ϕ = (ϕ_{1},…, ϕ_{I})^{T}, a popular version of the CAR model is characterized by

$$\varphi ~N(0,{[{D}_{\tau}(I-\alpha B)]}^{-1}),$$

(4)

where *B* is an *I* × *I* matrix with diagonal elements 0, and *D*_{τ} is a diagonal matrix with entries ${\tau}_{i}^{2}$, *i* = 1,…,*I*. A usual assumption is *D*_{τ} = τ^{2}*D*, where *D* is a diagonal matrix with common entries. On the basis of the formulation in equation (4), different choices of *B*, α (0, 1), and *D*_{τ} lead to various CAR modelling structures. Most often, we set *D* to be a diagonal matrix with entries *m _{i}, i* = 1,…,

For assessing practical benefits in modelling these effects, we also fit i.i.d models where ϕ ~ *N*(0, τ^{−2}*I*). Two basic schemes, first and last activation, of the cure rate model class should apply in both structures. We will analyse the resulting models and investigate the performances using Deviance Information Criterion (DIC; Spiegelhalter *et al.*^{33}). See section 4 for further details.

For investigating spatial associations in the latent link and the cure fraction jointly, the basic structure is similar to the univariate setting:

$$\begin{array}{r}\hfill \text{Model}\phantom{\rule{thinmathspace}{0ex}}1(\mathrm{a}):{\eta}_{ij}={\mathbf{x}}_{ij}^{\mathrm{T}}\beta +{\varphi}_{1i};\phantom{\rule{thinmathspace}{0ex}}\text{log}({\theta}_{ij})={\theta}_{0}+{\theta}_{2i};\phantom{\rule{thinmathspace}{0ex}}N~\text{Po}(\theta );\\ \text{Model}\phantom{\rule{thinmathspace}{0ex}}2:\text{log}({\theta}_{ij})\phantom{\rule{thinmathspace}{0ex}}={\mathbf{x}}_{ij}^{\mathrm{T}}\beta +{\varphi}_{2i};\phantom{\rule{thinmathspace}{0ex}}{\eta}_{ij}={\eta}_{0}+{\varphi}_{1i};\phantom{\rule{thinmathspace}{0ex}}N~\text{Po}(\theta ).\hfill \end{array}$$

Again, both schemes proposed by Cooner *et al.*^{16} should be identifiable. When multiple parameter sets (such as frailties in cure fractions and the Weibull link) need to be modelled jointly, as opposed to independently, we resort to multivariate CAR models originally proposed by Mardia^{34}; see also the works of Carlin and Banerjee^{35} and Gelfand and Vounatsou.^{36} Let ${\mathit{\nu}}_{i}^{\prime}=({\varphi}_{i1},\dots ,{\varphi}_{ip})$ be a vector of *p* variables associated with the *i*th region. Collecting these effects into $\mathit{\nu}\prime =({\mathit{\nu}}_{1}^{\prime},\dots ,{\mathit{\nu}}_{n}^{\prime})$, the joint distribution can be written down as

$$\mathit{\nu}~N(0,{[\mathbf{\Gamma}(I-{B}_{R})]}^{-1}),$$

(5)

where *B _{R}* is an

Mardia^{34} and Carlin and Banerjee^{35} discuss more general conditions for the positive-definiteness of **Γ**(*I* − *B _{R}*), although computational feasibility often encourages simpler specifications such as

Although the Kronecker structure offers computational and interpretational simplicity, there has been much research on developing more general spatial covariances, notably by Kim *et al.*,^{37} Carlin and Banerjee^{35} and Jin *et al.*^{20} The last work offers a Generalized Multivariate Conditionally Auto Regressive (GMCAR) model with emphasis on computational simplicity. Considering *p* = 2 and writing $\varphi \prime =({\varphi}_{1}^{\prime},{\varphi}_{2}^{\prime})$, where ${\varphi}_{1}^{\prime}=({\varphi}_{11},\dots ,{\varphi}_{n1})$ and ${\varphi}_{2}^{\prime}=({\varphi}_{12},\dots ,{\varphi}_{n2})$, the GMCAR models the joint distribution of ϕ as *p*(ϕ) = *p*(ϕ_{1}|ϕ_{2})*p*(ϕ_{2}). Writing this joint distribution as

$$\left(\begin{array}{c}{\varphi}_{1}\\ {\varphi}_{2}\end{array}\right)~N\left(\left(\begin{array}{c}0\\ 0\end{array}\right),\left(\begin{array}{cc}{\sum}_{11}& {\sum}_{12}\\ {\sum}_{12}^{\prime}& {\sum}_{22}\end{array}\right)\right),$$

(6)

setting $A={\sum}_{12}{\sum}_{22}^{-1}$ and ${\sum}_{11.2}={\sum}_{11}-{\sum}_{12}{\sum}_{22}^{-1}{\sum}_{12}^{\prime}$, we have ϕ_{1}|ϕ_{2} ~ *N*(*A*ϕ_{2}, Σ_{11.2}) and ϕ_{2} ~ *N*(0, Σ_{22}).

Treating ϕ_{1}|ϕ_{2} and ϕ_{2} as two univariate CARs, we assume ϕ_{1}|ϕ_{2} ~ *N*(*A*ϕ_{2}, τ_{1}[(*D* − α_{1}*W*)]^{−1}), and the marginal distribution for ϕ_{2} is *N*(0, τ_{2}[(*D* − α_{2}|*W*)]^{−1}). Furthermore, let *D* = Diag(*m _{i}*), and 0 < α

Our observed data comprises *D _{ij}* = (

$$L({D}_{ij};{\mathrm{\Omega}}_{ij},{N}_{ij})=P{(T\ge {t}_{ij}|{N}_{ij})}^{1-{\delta}_{ij}}\times {\left(-\frac{\mathrm{d}}{\mathrm{d}{t}_{ij}}P(T\ge {t}_{ij}|{N}_{ij})\right)}^{{\delta}_{ij}},$$

where *P*(*T* ≥ *t|N*) is as in equation (1), and $-(\mathrm{d}/\mathrm{d}t)P(T\ge t|N)=N\left(\begin{array}{c}N-1\\ r-1\end{array}\right){[S(t)]}^{N-r}$ [*F*(*t*)]^{r−1}*f* (*t*) with *r* = 1 or *N* accordingly for first or last activation. More generally, we can envision each subject-specific activation types, so that some of the *r*’s will be 1 and others will be *N*. However, it is more appropriate to think of the population (patients from a specific cancer type) as being subjected to a mechanistic generation of failure that determines the activation scheme for the first or last patient. So, we usually fix *r* = 1 or set it equal to *N* for each subject.

We seek the posterior distribution of **Ω** after marginalizing over the distribution of *N _{ij}*’s that considerably simplifies the MCMC implementation

$$L(\{{D}_{ij}\};\{{\mathrm{\Omega}}_{ij}\})={\displaystyle \prod _{i=1}^{I}{[{h}^{*}({t}_{ij})]}^{{\delta}_{ij}}{S}^{*}({t}_{ij}).}$$

(7)

Thus, we may sample *P*(**Ω**|{*D _{ij}*}) using say Metropolis–Hastings algorithms, without worrying about sampling the

Further computational simplifications arise as *h**(*t*) can often be evaluated in closed forms. For instance, in the first- and last-activation schemes. Under the first-activation scheme, we obtain *h**(*t*) = *h*(*t*)*E _{N}*[

For the last-activation scheme, we obtain the hazard as *h**(*t*) = *f* (*t*)*E _{N}*[

Finally, we offer model comparisons using the DIC^{33} as a measure of model choice. The DIC has nice properties for Gaussian likelihoods (as ours) and is particularly convenient to compute from posterior samples. This criteria is the sum of the Bayesian deviance (a measure of model fit) and the (effective) number of parameters (a penalty for model complexity). It rewards better fitting models through the first term and penalizes more complex models through the second term, with lower values indicating favourable models for the data. The deviance, up to an additive quantity not depending upon **Ω**, is simply minus twice the log-likelihood, *D*(**Ω**) = −2 log *L*({*D _{ij}*};

The National Cancer Institute’s SEER database, available online at http://seer.cancer.gov, provides a national cohort of women who have been monitored for assessing breast cancer prognosis. In addition, available individual-level covariates are age at diagnosis, the number of primary cancers each woman has had diagnosed, the stage of the disease (local (200), regional (94 patients) or distant (4 patients), with local as baseline), and the county information at the time of diagnosis. We consider a sample of 298 patients from Iowa (with *99* counties), who were all diagnosed of breast cancer in January 1992 and monitored through 1998. The response here is time to death from breast cancer: only those who have been identified as having died from metastasis of nodes in the breast (there were 49 such deaths) are considered having failed, while the rest (including those who might have died from other types of cancer or other causes) are considered censored. The longest observation time of this sample is 83 months.

Under both activation schemes in our analysis of the cancer data, we assumed that the regression coefficients β had a non-informative *N*(0, 10^{3}) prior (flat priors are admissible here as well), while a relatively weak Gamma(2, 0.001) prior (mean 2000) was used for ρ. Specifically, for model 2, an *N*(0, 100) prior for η was used, a Gamma(2, 0.001) prior was used for θ in model 1(a), while *U*(0, 1) priors were assigned to θ in models 1(b) and (c). We also use Gamma(2, 0.001) for the precision τ in either independent case or the CAR model. These priors are vague enough to allow the data to drive the posterior inference, while still leading to acceptable MCMC convergence. We experimented with other hyperparameter values and did not observe much sensitivity to these choices in the posterior distributions.

For each analysis, we ran two initially dispersed parallel MCMC chains for 30 000 iterations each, where convergence was monitored using sample autocorrelations within the chains, cross-correlations between the parameters, and plots of the sample traces. These tools suggested discarding the first 5 000 iterations from each chain as pre-convergence burn-in. Retaining every 10th of the remaining 2 × 25 000 = 50 000 iterations yielded a final sample of size 5 000 for posterior analysis. The first- and last-activation models were implemented in `WinBUGS` (see www.mrc-bsu.cam.ac.uk/bugs/welcome.shtml) using R (www.r-project.org) for final posterior summarization. The relevant codes are available from www.biostat.umn.edu/~sudiptob/Software/Software.html.

We use univariate spatial models to analyse the spatial patterns in the latent link η and the cure fraction parameter θ separately. We use a single set of frailties $({\left\{{\varphi}_{1}\right\}}_{i=1}^{I})$ and fit models 1(a)–(c) regressing in η and model 1(a) and model 2 with regression in log(θ). Table 1 presents the DIC scores for these models when the ϕ_{i}’s are i.i.d. Gaussian random effects. For regression in η, we see that model 1(c) has the best DIC scores in both the first- and last-activation models. Interestingly, first-activation scores show slightly greater variation between the models in their DIC scores as compared to the last-activation model. This might reflect greater sensitivity of the first activation to the distributional assumptions of *N*. Also, for model 1(c) the effective number of parameters, *p _{D}*, seem to be prominently higher due to greater variation in the random effects. For regression in log(θ), we find that model 2 has a substantially lower DIC score compared to model 1(a) under the first-activation scheme, while the latter performs moderately better than the former for the last activation setting. A very similar story is obtained from Table 2, where the i.i.d. frailties are replaced by spatially associated CAR frailties.

Univariate analysis: DIC values for various models under first- and last-activation schemes and frailties (with i.i.d. normal assumption) in latent link η or cure fraction parameter θ. The model with the best DIC value in each column is **...**

Univariate analysis: DIC values for various models under first- and last-activation schemes and spatial frailties in latent link η or cure fraction parameter θ with CAR assumption

We did not observe much significance for the covariates included in our model. Indeed, the incorporation of spatial effects sometimes widen these credible intervals, as compared to models without random effects. Our preliminary results apparently indicate a rather slim impact of the activation schemes on the posterior estimation of spatial clustering. Therefore, we present the maps of the spatial posterior medians for η and θ only for model 1(a) and model 2 under first-activation schemes. Also, for η the different structures in model 1 did not cause major differences in spatial patterns, while for θ only model 1(a) and model 2 yielded stable convergence – a fact corroborated by the theoretical impropriety of the others. More specifically, the i.i.d. frailties, plotted in Figures 1 and and2,2, show less pronounced spatial patterns. Both figures reveal only slight differences in spatial patterns for the cure rates or the latent mean between model 1(a) and model 2. Also, we can see that the large difference of the ϕ values in model 1(a) and model 2, suggesting that the simple model (univariate and i.i.d. normal on ϕ) are not stable or reliable (perhaps due to the small sample sizes in some counties).

Univariate analysis – i.i.d. normal (frailties in latent link η): maps for median frailties under first-activation scheme for model 1(a) (left) and model 2 (right).

Univariate analysis – i.i.d. normal (frailties in cure fraction parameter θ): maps for median frailties under first-activation scheme for model 1(a) (left) and model 2 (right).

Figures 3 and and44 reveal much clearer spatial patterns from the northeast to the southwest for both η and θ. Since the latent mean survival time and the cure fractions are non-increasing functions of η and θ respectively, the maps for all the models suggest that the southwest has a better survival experience than the northeast. Although both the southwest and the northeast are relatively less populated rural regions in Iowa, these clusters might reflect unknown or missing covariates that explain differences in patient care (e.g., access to health care, quality of medical facilities and so on). Incorporating regressors containing such information might well lead to better smoothing of the maps.

Univariate analysis – CAR (spatial frailties in latent link η): maps for median spatial frailties under first-activation scheme for model 1(a) (left) and model 2 (right).

Here we consider joint spatial associations in η and θ simultaneously. In Table 3 we show the different models we consider here and their DIC scores. We classify the models as those having two sets of i.i.d. random effects, those with two independent CAR distributions and finally with jointly associated GMCAR random effects. Not much variation is seen across the DIC scores. Under the first-activation scheme we find model 2 to be slightly better for the i.i.d. effects and the two independent CAR effects, while model 1 seems to be marginally better for the GMCAR effects. The story is somewhat similar for the last-activation setting, although we did not find suitable convergence for model 2 from the data under the GMCAR setting.

Multivariate analysis: DIC values for various models under first- and last-activation schemes, with i.i.d. normal, CAR or GMCAR assumption on spatial frailties

In Tables 4–6, we present parameter inference results from the best model under first and last activation, with three bivariate assumptions on spatial frailties. For all six models, we find insignificant impact for age at diagnosis and the number of primary nodes, although age at diagnosis has a positive median for each of the models. As expected, the stage variable, indicating the extent of the disease, “Regional” and “Distant” lead to significant hazard increase relative to “Local”. Finally, the Weibull scale parameter ρ is robustly estimated across the models.

Posterior quantiles for different models (corresponding to best DIC) under first- and last-activation schemes, with i.i.d. normal assumption on spatial frailties

Posterior quantiles for different models (corresponding to best DIC) under first- and last-activation schemes, with GMCAR assumption on spatial frailties

Turning to the spatial maps, two sets of i.i.d. Gaussian random effects in η and θ do not bring about much changes in the posterior behaviour for ϕ, so we do not present the maps for space. However, in Figures 5 and and6,6, considering the spatial frailties in the frailties together, we find that spatial clustering with frailties in η is inconsistent with that in log(θ). This happens for model 1(a) and model 2 and is in clear contrast with Figures 3 and and44 which displayed consistent behaviour. The differences between two models for both parameters are large and the spatial pattern seem reversed. It is worth mentioning that the values now are very well stabilized. Two separate i.i.d. normal or CAR frailties might still not be enough to capture the spatial patterns.

Bivariate analysis – two independent CAR (spatial frailties in latent link η i.e., ϕ_{1}): maps for median spatial frailties under first-activation scheme for model 1(a) (left) and model 2 (right).

Bivariate analysis – two independent CAR (spatial frailties in cure fraction parameter ϕ i.e., ϕ_{2}): maps for median spatial frailties under first-activation scheme for model 1(a) (left) and model 2 (right).

These lead us to question some of the fundamental assumptions for the independence between the spatial effects in the latent link η and the cure fraction parameter θ. Simply introducing two independent univariate spatial structures is inadequate, hence the MCAR or the GMCAR should be implemented into our models. With the GMCAR assumption, the maps (Figures 7 and and8)8) show somewhat similar spatial trends as the univariate cases, which may suggest a lack of significant spatial association for this set of data.

Bivariate analysis – GMCAR (spatial frailties in latent link η i.e., ϕ_{1}): maps for median spatial frailties under first-activation scheme for model 1(a) (left) and model 2 (right).

We plan to extend our current work into several directions. First, we will substantially enhance the current illustration with the hierarchical activation scheme results. Next we propose to implement these models with more flexible MCAR specifications. Finally, we will investigate spatio-temporal effects to account for temporal associations where we may model the cure fractions θ_{i}(*t*) or the latent link η_{i}(*t*) as functions of time. These may include direct time-varying regression, polynomial or spline functions of time, or even completely non-parametric modelling using temporal processes.

Freda Cooner, Division of Biostatistics, School of Public Health, University of Minnesota, Minneapolis, MN, USA.

Sudipto Banerjee, Division of Biostatistics, School of Public Health, University of Minnesota, Minneapolis, MN, USA.

A Marshall McBean, Division of Health Services Research and Policy, School of Public Health, University of Minnesota, Minneapolis, MN, USA.

1. Lawson AB, Biggeri A, Boehning D, Lesaffre E, Viel J-F, Bertollini R. Disease mapping and risk assessment for public health. Wiley; 1999.

2. Elliot P, Wakefield JC, Best NG, Briggs DJ. Spatial epidemiology: methods and applications. Oxford University Press; 2001.

3. Lawson AB, Williams FLR. An introductory guide to disease mapping. Wiley Medical Sciences; 2001.

4. Banejee S, Wall MM, Carlin BP. Fraity modelling for spatially correlated survival data with application to infant mortality in Minnesota. Biostatistics. 2003;4:123–142. [PubMed]

5. Li Y, Ryan L. Modeling spatial survival data using semiparametric frailty models. Biometrics. 2002;58:287–297. [PubMed]

6. Cox DR, Oakes D. Analysis of survival data. Chapman & Hall; 1984.

7. Berkson J, Gage RP. Survival curve for cancer patients following treatment. Journal of the American Statistical Association. 1952;47:501–515.

8. Farewell VT. The use of mixture models for the analysis of survival data with long term survivors. Biometrics. 1982;38:1041–1046. [PubMed]

9. Farewell VT. Mixture models in survival analysis: are they worth the risk? Canadian Journal of Statistics. 1986;14:257–262.

10. Goldman AI. Survivorship analysis when cure is a possibility: a Monte Carlo study. Statistics in Medicine. 1984;3:153–163. [PubMed]

11. Ewell M, Ibrahim JG. The large sample distribution of the weighted log-rank statistic under general local alternatives. Lifetime Data Analysis. 1997;3:5–12. [PubMed]

12. Yakovlev AY, Asselain B, Bardou VJ, Fourquet A, Hoang T, Rochefordiere A, Tsodikov AD. A simple stochastic model of tumor recurrence and its application to data on premenopausal breast cancer. In: Asselain B, Boniface M, Duby C, Lopez C, Masson JP, Tranchefort J, editors. Biometrie et Analyse de Donnees Spatio-Temporelles, No. 12. Société Française de Biométrie; 1993. pp. 66–82.

13. Chen M-H, Ibrahim JG, Sinha D. A new Bayesian model for survival data with a surviving fraction. Journal of the American Statistical Association. 1999;94:909–919.

14. Ibrahim JG, Chen M-H, Sinha D. Bayesian survival analysis. Springer-Verlag; 2001.

15. Tsodikov A, Ibrahim J, Yakovlev A. Estimating cure rates from survival data: an alternative to two-component mixture models. Journal of the American Statistical Association. 2003;98:1063–1078. [PMC free article] [PubMed]

16. Cooner F, Banerjee S, Carlin BP, Sinha D. Research Report 2005–06. Division of Biostatistics, University of Minnesota; Flexible cure rate modelling under latent activation schemes.

17. Banerjee S, Carlin BP. Parametric spatial cure rate models for interval-censored time-to-relapse data. Biometrics. 2004;60:268–275. [PubMed]

18. Carlin BP, Louis TA. Bayes and empirical Bayes methods for data analysis. second edition. Chapman & Hall/CRC Press; 2000.

19. Rue H, Held L. Gaussian Markov Random Fields; Theory and Application. Chapman & Hall/CRC Press; 2005. Monographs in Statistics and Applied Probability No. 104.

20. Jin X, Carlin BP, Banerjee S. Generalized hierarchical multivariate CAR models for areal data. Biometrics. 2006 in press. [PubMed]

21. Rao CR. Linear statistical inference and its applications. second edition. Wiley; 1973.

22. Tucker SL, Thames HD, Taylor JMG. How well is the probability of tumor cure after fractionated irradiation described by Poisson statistics? Radiation Research. 1990;124:273–282. [PubMed]

23. Tucker SL, Taylor JMG. Improved models of tumor cure. International Journal of Radiational Biology. 1996;70:539–553. [PubMed]

24. Moolgavkar SH, Luebeck EG, De Gunst MC. Two mutation model for carcinogenesis: relative roles of somatic mutations and cell proliferation in determining risk. In: Moolgavkar SH, editor. Scientific issues in quantitative cancer risk assessment. Birkhauser; 1990. pp. 136–152.

25. Besag J. Spatial interaction and the statistical analysis of lattice system (with discussion) Journal of the Royal Statistical Society, Series B. 1974;36:192–236.

26. Cressie NAC. Statistics for spatial data. revised edition. Wiley; 1993.

27. Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data. Chapman & Hall/CRC Press; 2004.

28. Heikkinen J, Högmander H. Fully Bayesian approach to image restoration with an application in biogeography. Applied Statistics. 1994;43:569–582.

29. Högmander H, Møller J. Estimating distribution maps from atlas data using methods of statistical image analysis. Biometrics. 1995;51:393–404.

30. Hoeting JA, Leecaster M, Bowden D. An improved model for spatially correlated binary response. Journal of Agricultural, Biological, and Environmental Statistics. 2000;5:102–114.

31. Besag J, York JC, Mollié A. Bayesian image restoration, with two applications in spatial statistics (with discussion) Annals of the Institute of Statistical Mathematics. 1991;43:1–59.

32. Lawson AB, Browne WJ, Vidal Rodeiro CL. Disease mapping with WinBUGS and MLwiN. Wiley; 2003.

33. Spiegelhalter DJ, Best N, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit (with discussion) Journal of the Royal Statistical Society, Series B. 2002;64:583–639.

34. Mardia KV. Multi-dimensional multivariate Gaussian Markov random fields with application to image processing. Journal of Multivariate Analysis. 1998;24:265–284.

35. Carlin BP, Banerjee S. Hierarchical multivariate CAR models for spatio-temporally correlated survival data (with discussion) In: Bernardo JM, Bayarri MJ, Berger JO, Dawid AP, Heckerman D, Smith AFM, West M, editors. Bayesian Statistics. Oxford University Press; 2003. pp. 45–63.

36. Gelfand AE, Vounatsou P. Proper multivariate conditional autoregressive models for spatial data analysis. Biostatistics. 2003;4:11–25. [PubMed]

37. Kim H, Sun D, Tsutakawa RK. A bivariate Bayes method for improving the estimates of mortality rates with a twofold conditional autoregressive model. Journal of the American Statistical Association. 2001;96:1506–1521.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |