TEXT SIZE

CrossRef (0)
Cure rate proportional odds models with spatial frailties for interval-censored data

Correspondence to: 1Corresponding author: Departamento de Matemática Aplicada e Estatística, Universidade de São Paulo, Avenida Trabalhador São-carlense, 400 - Centro CEP: 13566-59, São Carlos-SP, Brasil. E-mail: suzuki@icmc.usp.br
Received July 8, 2017; Revised September 26, 2017; Accepted November 7, 2017.
Abstract

This paper presents proportional odds cure models to allow spatial correlations by including spatial frailty in the interval censored data setting. Parametric cure rate models with independent and dependent spatial frailties are proposed and compared. Our approach enables different underlying activation mechanisms that lead to the event of interest; in addition, the number of competing causes which may be responsible for the occurrence of the event of interest follows a Geometric distribution. Markov chain Monte Carlo method is used in a Bayesian framework for inferential purposes. For model comparison some Bayesian criteria were used. An influence diagnostic analysis was conducted to detect possible influential or extreme observations that may cause distortions on the results of the analysis. Finally, the proposed models are applied for the analysis of a real data set on smoking cessation. The results of the application show that the parametric cure model with frailties under the first activation scheme has better findings.

Keywords : Bayesian inference, influence diagnostics, survival models with a cure fraction, spatial frailty
1. Introduction

Advances in medical and health science have led to several clinical studies that showed non-negligible proportions of patients who respond favorably to a treatment and are non-susceptible to the event of interest. Such a proportion of non-susceptible patients is considered cured or with prolonged disease-free survival which is usually referred to as the cured fraction. Survival models with cure fraction, also known as long-term survival models or cure rate models, have been widely developed to accommodate a cured fraction.

One of the most popular cure rate model is the mixture model introduced by Berkson and Gage (1952). It assumes that a certain proportion of the patients are cured, in the sense that they do not present the event of interest during a long period of time and can be seen as cured or immune from the causes of death under study (Rodrigues et al., 2009). Mixture models can also accommodate the effects of covariances on the cure rate of populations and the survival function of susceptible to (non-cured) individuals via two components of separate regression. An alternative was later proposed by Yakovlev and Tsodikov (1996) and Chen et al. (1999) concerning cancer therapy studies. In this context, it is assumed that, for each individual, there is a latent quantity M of cells able to cause cancer, being also assumed that M has a Poisson distribution. The model is known as the promotion time cure model or bounded cumulative hazard model, where an inference can be made concerning the distribution of cells that produce cancer. This subject has attracted the attention of several researchers, because it adheres to the proportional hazard structure and incorporates parameters with clear biological interpretations.

Gu et al. (2011) have recently developed another class of cure rate models where the survival function has a proportional odds structure that is an already popular regression model in ordinary survival model (Bennett, 1983; Collett, 1994). The proportional odds models with cure rate retain all the advantages of regular proportional odds models regarding the survival analysis. They can also be derived from the latent factors models conducted by Cooner et al. (2007), considering these models share some properties.

In many clinical trial studies, it is very common for patients to be periodically examined for disease occurrence or progression. In such a situation, the exact failure time which patients cannot be observed exists, but can only be allowed to lie in an interval after a sequence of examination times that is known as interval-censored (Peto, 1973). The estimation methods available to the right-censored data, such as the Kaplan-Meier estimator, is not considered adequate to be applied in interval-censored data, because it can lead to biased estimation and invalid inferences. The information concerning interval censorship should be taken into account in modeling (Lindsey and Ryan, 1998; Rücker and Messerer, 1988; Sun and Chen, 2010).

The lifetimes data are sometimes collected from several regions, which can lead to different effects for each observation. To consider these different effects in the cure model for interval-censored data, Xiang et al. (2011) adopted the mixture model by including two random effects (also known as frailties) to each cluster to explain the effects on the survival time for susceptible individuals and the effects on the cure fraction. Hu and Xiang (2013) adopted the non-mixture cure model and introduced random effects in the promotion of the cure time model. The development of the geographic information system has allowed several researchers to incorporate geographical information on the subjects under study. They have also developed survival models that account for spatial clustering and variation. Banerjee and Carlin (2004) have specified spatial correlation to fit geographically clustered interval-censored survival data using a parametric cure rate model. They have developed a Bayesian approach to the mixture cure model with spatial random effects on the survival function for subjects at risk and spatial frailties by using a multivariate conditionally autoregressive (MCAR) prior. Pan et al. (2014) propose a Bayesian approach under a proportional hazards frailty model to analyze interval-censored survival data with spatial correlation.

This paper analyzes smoking cessation data. In a smoking cessation study, all patients (smokers) were randomized both into a smoking intervention (SI) group, or into a usual care (UC) group that receives no special anti-SI. The SI program treatments were conducted in Rochester city, being located in the center of the map. The details concerning the programs can be found in Murray et al. (1998). Here, each patient was observed once a year over the period of five years of monitoring. Our event of interest concerns if they will relapse into smoking (resume smoking) or not. If a smoker starts smoking again after an initial attempt to quit, then only an approximate one-year time interval was observed from the previous observation to the current observation. Thus, the relapse times are interval-censored. In this analysis, we limit our attention to patients who are known to have quit smoking at least once during the period of study and who have an identifiable Minnesota zip code of residence. Thus, the data consists of a total of 223 patients who reside in 51 zip codes in the Southeastern corner of Minnesota; among them, there are 65 patients that have underwent a relapse, which implies the empirical cure rate is approximately 71%. Ma and Xiang (2013) also confirmed the existence of a non-eligible cure fraction in the population.

In this paper, we present a proportional odds cure model (Gu et al., 2011) to allow spatial correlations by including spatial frailty in the interval-censored data setting, to propose parametric cure rate models that include independent and dependent spatial frailties. These new cure rate models enable different underlying activation mechanisms which lead to the event of interest. Also, we develop inference procedures through a Bayesian perspective. The proposed models can be also compared to models introduced by Banerjee and Carlin (2004) and Pan et al. (2014) through the deviance information criterion (DIC) proposed by Spiegelhalter et al. (2002).

We also conduct influence diagnostics to examine the assumptions of the model and conduct studies on sensitivity to detect possible influential or extreme observations that can cause distortions in the results of the analysis. Here the influence of the diagnostics of case deletion are developed for the posterior joint distribution based on the ψ-divergence (Peng and Dey, 1995; Weiss, 1996).

The remainder of the paper is organized as follows. In Section 2, the introduction of frailty models and some distributions for spatial frailties are presented and followed by the spatial distributional aspects of our modeling. The Bayesian inference for the proposed models is developed in the Section 3. In Section 4, the proposed models are fit into a read data set (smoking cessation study). Finally, Section 5 concludes with some general remarks.

2. Cure rate proportional odds model

Gu et al. (2011) have recently proposed a cure rate proportional odds model (CRPO model). Unlike the model proposed by Chen et al. (1999), the ratio of hazards for the CRPO model for two covariate values does not remain over time. The model can be characterized by the latent factors model by Cooner et al. (2007), which contains a Geometric distribution for the number of latent factors. Supposing that there are I regions and ni individuals in ith region. Let Ti j denotes the random variable for time to the event of the jth individual in the ith region, where j = 1, …, ni and i = 1, …, I. We suppose that the (i, j)th individual is potentially exposed to Mi j latent risks, where Mi j denote the initial number of competing causes concerning the occurrence of an event, and assuming Mi j has a Geometric distribution with parameter 1/(1 + θi j), the probability mass function is given by

$P(Mij=m)=θijm(θij+1)m+1, m=0,1,2,…,$

where θi j > 0, E(Mi j) = θi j, and Var(Mi j) = θi j(1 + θi j).

Let Yci j denotes the lifetime (for relapse into smoking) of jth individual in ith region due to the cth (c = 1, …, Mi j) latent cause or risk. Common latent causes of relapse into smoking are giving in to cravings, stress, weight gain and nicotine withdrawal symptoms, which may include dizziness, nausea, constipation, increased appetite, anxiety, depression, and sensation of tightness in the chest (National Cancer Institute at http://smokefree.gov/withdrawal). The problem is that the causes or risks are latent in the sense that it is hard to know the exact cause of relapse into smoking.

Given Mi j = m, Y1i j, Y2i j, …, Ymi j are assumed to be independent and identically distributed with a common distribution function F(·) = 1 − S (·) that does not depend upon Mi j. If we assume that the presence of any latent risk will ultimately lead to the occurrence of the event, the time to the event of interest Ti j could be defined as Ti j = min{Y1i j, …, YMi ji j} for Mi j ≥ 1. In this case we are assuming the cause, which comes first, not even knowing what it is, cause relapse into smoking. If Mi j = 0, then the individual is not at risk of occurrence of the event (relapse into smoking)) and is considered cured. In this case, we define Ti j = ∞ with P(Ti j = ∞|Mi j = 0) = 1. Thus, the survival function for the population is given by

$Spop(tij)=[1+θijF (tij)]-1.$

Note that this survival function has a proportional odds structure when covariates xi j are modeled via θi j(xi j) and the latent survival F(ti j) is free of xi j, because

$1-Spop(tij∣xij)Spop(tij∣xij)=θij(xij)F (tij).$

This CRPO model can also be obtained as an extension of the transformation-cure model (Zeng et al., 2006).

The probability density function (pdf) and the hazard function associated to (2.2) are given by

$fpop(tij)=θijf(tij)[1+θijF (tij)]-2 and hpop(tij)=θijf(tij)[1+θijF (tij)]-1,$

respectively, where f (ti j) = (∂/∂ti j)F(ti j).

Note that, the survival function in (2.2) can also be written as a mixture cure model

$Spop(tij)=(1+θij)-1+(1-(1+θij)-1)[{1+θijF (tij)}-1-(1+θij)-11-(1+θij)-1].$

Thus, the survival functions of uncured (susceptible) individuals can be expressed by

$Ssus(tij)=[1+θijF (tij)]-1-(1+θij)-11-(1+θij)-1.$

Now, if we assume another situation where the presence of all latent risks will ultimately lead to the occurrence of the event. In this case we are assuming the cause, which is later, not even knowing what it is, cause relapse into smoking. Thus, the time to the event of interest is defined by the random variable Ti j = max{Yci j, c = 1, …, Mi j} for Mi j ≥ 1 and Ti j = ∞ if Mi j = 0 with P(Ti j = ∞|Mi j = 0) = 1. The survival function for the population is given by

$Spop(tij)=1+(1+θij)-1-[1+θijS(tij)]-1.$

The corresponding pdf and the hazard function are given by

$fpop(tij)=θijf(tij)(1+θijS(tij))-2,$

and

$hpop(tij)=θf(tij)[1+θijS (tij)]-21+(1+θij)-1-(1+θijS (tij))-1,$

respectively. The survival function (2.3) can also be written as a mixture cure model

$Spop(tij)=(1+θij)-1+(1-(1+θij)-1)1-(1+θijS (tij))-11-(1+θij)-1.$

Thus, the survival functions of susceptible individuals is given by

$Ssus(tij)=1-(1+θijS (tij))-11-(1+θij)-1.$

The first situation is also known as first activation (FA) scheme because, in this case, we assume the event of interest occurs when the first possible cause is activated. However, the second situation is known as the last activation (LA) scheme because the event of interest only takes place after all the latent causes have been activated (Cooner et al., 2007). Thus, we denoted the survival functions (2.2) and (2.3) by $SpopF(tij)$ and $SpopL(tij)$, respectively. There is another kind of situation where the event of interest occurs: when some of the possible causes are activated and, given the number of latent causes Mi j, the number of activated causes is a random variable with discrete Uniform distribution {1, …, Mi j}. This situation is known as random activation scheme. In this case, the survival function for the population is given by

$SpopR(tij)=(1+θij)-1+(1-(1+θij)-1)S (tij),$

where the superscript R denotes random activation scheme.

Note that whichever the activation scheme, the density and hazard functions of the cure models are improper functions, since the survival functions are not proper. Its cure fraction is the same for all the activation schemes and can be obtained by p0i j = limti j→∞Spop(ti j) = (1 + θi j)−1. However, under different activation schemes, the models differ by its survival, density, and hazard functions. Moreover, under the conditions of the models (2.2), (2.3), and (2.4) for any distribution function F(·), we have $SpopF(tij)≤SpopR(tij)≤SpopL(tij)$ for all ti j > 0.

The cure fraction plays a key role in the survival models with cure fraction. So we consider the parametrization of the model in terms of the cure fraction. Since p0i j = (1 + θi j)−1, we have $θij=p0ij-1-1$. Moreover, we propose that the cured fraction of an individual (i, j)th be associated with covariates xi j. Thus linking p0i j to covariates xi j by

$p0ij=exp(ξij)1+exp(ξij), j=1,…,ni, i=1,…,I,$

where ξi j is a linear function of covariates, $ξij=xij′b$ with b is a p1-dimensional vector which represents the effects of covariates on the cured fraction. Thus, the models (2.2), (2.3), and (2.4) parametrized in terms of p0i j can be written as

$SpopF(tij)=[1+(p0ij-1-1)F (tij)]-1,$$SpopL(tij)=1+p0ij-[1+(p0ij-1-1)S (tij)]-1$$SpopR(tij)=p0ij+(1-p0ij)S (tij).$

The model in (2.7) is the same considered by Banerjee and Carlin (2004).

Henceforth, we assume the Weibull distribution for the latent variables Yci j, c = 1, …, i = 1, …, I, and j = 1, …, ni in (2.2), (2.3), and (2.4) with $F(ycij∣φ)=exp(-ycijαeλij)$ where φ = (α, λi j), α > 0 is a shape parameter and λi j ∈ ℝ is a scale parameter with $λij=zij′β$ is the linear predictor of the covariates, where zi j represents the covariates of an individual (i, j) and β is the p2-dimensional vector that represents the effects of covariates on the survival model component. This parametrization of the Weibull distribution is implicitly assumed as a proportional hazards model with the baseline hazard function h0(t|α) = αtα−1. Here, we called the model in (2.2) the Weibull cure rate model (WCRM) under the FA denoted by WCRM-FA and (2.3) by WCRM under the LA denoted by WCRM-LA.

By following Banerjee and Carlin (2004), we introduce the frailties Ui and Vi to evaluate the effect of regions on the lifetime of individuals and on the cured fraction through linear predictor

$λij=zij′β+Ui,ξij=xij′b+Vi, for j=1…,ni, i=1,…,I.$

Here, the frailties Ui and Vi are spatially correlated across the regions. In this work, we propose two approaches. Firstly, we employ separate independent conditionally auto-regressive (CAR) prior distribution on (U, V).

The CAR model was originally developed by Besag (1974). These priors have become very popular in the Bayesian analysis of real data, especially in disease mapping. Let ψ1 = U and ψ2 = V denote the random frailties vectors, W denote the adjacent matrix of the map, so that Wii = 1 if the regions i and i′ are adjacent, otherwise it would be 0, while wi+ = ∑j wi j denotes the number of regions adjacent to region i. The zero-centered CAR prior for ψl, l = 1, 2 is an I-dimensional Gaussian distribution with mean 0 and precision matrix $ρl-1(DW-W)$, where ρl > 0 and DW is diagonal with (DW)ii = wi+ for i = 1, …, I. For l fixed, the CAR model is denoted by CAR(ρl). Note that the precision matrix $ρl-1(DW-W)$ is rank deficient and is therefore a non-positive definite matrix that leads the distribution to improper. This singularity, although theoretically awkward, it creates a small problem during the Bayesian implementation, since the identifying sum-to-zero constraint $∑i=1Iψil=0$ is easily imposed in a Gibbs sampler simply by re-centering the φi draws around zero after every iteration.

Second, we assume that the spatial priors on (U, V) are dependent, presenting MCAR prior distribution. Let ψ = (U′, V′)′, we employ a multivariate MCAR distribution with a common smoothness parameter a, that is, ψ has a normal distribution with mean 0 and precision matrix Λ ⊗ (DWaW), where ⊗ denotes the Kronecker product, Λ is a 2×2 symmetric and positive definite matrix, a ∈ (0, 1) and W is standardized so that each of its rows sum to 1. This prior is a proper distribution, the parameter a has a spatial smoothness interpretation. The value of a closer to 1 implies greater weight on the adjacency matrix W, while a close to 0 implies the adjacency structure has few roles to play in the precision matrix. This prior is denoted by MCAR(a, Λ).

We also assume the parameter ψ with an extended MCAR distribution proposed by Gelfand and Vounatsou (2003) and Carlin and Banerjee (2003), which assumes different smoothness parameters for the parameters U and V, namely a1 and a2. Let (DWaiW) denotes the corresponding positive definite matrix and $Ri′Ri$ denotes its Cholesky factorization, where Ri is an upper triangular matrix with real and positive diagonal entries with I × I dimension for i = 1, 2. The precision matrix can be written as

$Λ⊗(DW-aW)=[λ11R1′R1λ12R1′R2λ21R2′R1λ22R2′R2]=R′(Λ⊗II)R,$

where λi j are the elements of Λ. The $R=[R100R2]$ has dimension 2I × 2I, and the Λ ⊗ (DWaW) is positive definite since Λ is positive definite. We denoted this prior by MCAR(a1, a2, Λ).

3. Bayesian inference

Let = {(Ai j, xi j, zi j, δi j); j = 1, …, ni, i = 1, … M} denote the observed data, where Ai j = (ti jL, ti jR] is the interval during which individual j in cluster i occur the event of interest, xi j and zi j are the p1-dimensional and p2-dimensional vectors of covariates, and δi j is following interval censoring indicator: δi j = I(ti jR < ∞). For a special case where the survival time is right (left) censored, Ri j = +∞ (Li j = 0), whereas for exact observations, it is ti jL = ti jR. By following Finkelstein (1986), the likelihood function regarding the conditional distribution of given (U, V) for the general interval-censored cure rate model is given by

$L(ϕ∣D,U,V)∝∏i=1I∏j=1ni(Spop(tijL∣ϕ)-Spop(tijR∣ϕ))δij Spop(tijL∣ϕ)1-δij∝∏i=1I∏j=1niSpop(tijL∣ϕ)(1-Spop(tijR∣ϕ,)Spop(tijL∣ϕ))δij,$

where ϕ = (b, β, U, V, α) and α is the shape parameter of the Weibull distribution. For a Bayesian analysis, we assume the follow prior densities for parameters the b′, β′, and α

• $bj~N(μb,σb2)$, j = 0, …, (p1 − 1), with μb and σb known;

• $βj~N(μβ,σβ2)$, j = 1, …, p2, with μβ and σβ known;

• $α~N(μα,σα2)I(0,∞)$, with μα and σα known;

where N(μ,σ2)I(a,b) denotes the truncated normal distribution that is the probability distribution of a normally distributed random variable whose value lies within the interval −∞ ≤ a < b ≤ ∞. In several areas, especially in medicine, it is preferable to use the prior information when it is available. In addition, a truncated normal distribution as prior facilitates the insertion of information in certain regions of the parameter space, since the hyperparameters no longer represent the mean and variance, but still controls the region of higher probability mass.

### 3.1. Independent assumption

For the independent assumption, we employ the separate independent CAR prior on the random frailties U = (U1, …., UI )′ and V = (V1, …, VI )′, that is,

• U1, …, UI ~ CAR(ρ1);

• V1, …, VI ~ CAR(ρ2);

where ρ1 and ρ2 are positive unknown hyper-parameters. We assume they have Inverse-Gamma prior with the known shape parameter a0 > 0 and scale parameter b0 > 0. Thus, the joint posterior distribution for (ϕ, ρ1, ρ2) is given by

$π(ϕ,ρ1,ρ2∣D)∝L(ϕ∣D)π(ϕ,ρ1,ρ2),$

where π(ϕ, ρ1, ρ2) = π(b)π(β)π(α)π(U|ρ1)π(V|ρ2)π(ρ1)π(ρ2) is the joint prior distribution for (ϕ, ρ1, ρ2) and is the likelihood function given in (3.1). Note that, this joint posterior density is analytically intractable. Thus, we based our inference on Markov chain Monte Carlo (MCMC) simulation methods. We can observed that the full conditional distributions for the parameters b, β, α, U while V presents no closed forms; therefore, we will use the Metropolis-Hastings algorithm to generate posterior samples for these parameters. To avoid range restrictions concerning the parameter α, we define ζ = log(α) to transform all parameters space into real space (necessary to work with the densities from the Gaussian proposal). Let ϑ = (b, β, ζ, ρ1, ρ2), using the Jacobian of this transformation, the joint prior density π(ϑ) can be expressed by

$π(ϑ)=π(b,β,ζ-1,U,V,ρ1,ρ2)×exp(∑i=1Qζi),$

where ζ−1 denotes the inverse function of ζ, that is, $ζ-1={ζi-1=exp(ζi), i=1…,Q}$.

However, the full conditional distributions for parameters ρi are given by

$π(ρi∣ϑ-ρi,D)∝π (ψi∣ρi) π(ρi)∝(ρi)-k2exp(-12ρiψi′(DW-W)ψi)ρi-a0-1 exp (-b0ρi-1)∝ρi-(a0+k2)-1exp{-(ψi′(DW-W)ψi2+b0)ρi-1}, i=1,2,$

where ψ1 = U, ψ2 = V and k is the rank of the matrix DWW. Thus, the full conditional distributions of the parameter ρi is an Inverse-Gamma distribution with parameters $a0+(k/2) e b0+(1/2)(ψi′(DW-W)ψi)$. In this case, the Gibbs sampler algorithm (see Gamerman and Lopes, 2006) is used to generate a posterior sample.

Thus, the joint posterior density of π(ϑ|D) is proportional to

$L(ϕ)=exp{-12[σb-2∑i=0p1-1bi2+σβ-2∑i=1p2βi2+∑i=1Qexp(2ζi)σα2+U′(DW-W)Uρ1+V′(DW-W)Vρ2]-(a0+1)(log(ρ1)+log(ρ2))-(b0ρ1+b0ρ2)+∑i=1Qζi}.$

### 3.2. Dependent assumption

Now we assume the spatial priors of the parameters (U, V) dependent from each other. Let ψ = (U′, V′)′, we first employed the parameter ψ as a MCAR distribution with a common smoothness parameter a, i.e.,

$ψ~MCAR(a,Λ).$

Further, we employ the parameter ψ an extend MCAR distribution that assumes the different smoothness parameters for the parameters U and V, namely a1 and a2, that is,

$ψ~MCAR(a1,a2,Λ).$

Here, we consider the prior distributions for the parameters a and Λ the same that were used by Carlin and Banerjee (2003), that is,

• ai ~ Uniform(0, 1) or ai ~ Beta(18, 2), for i;

• Λ ~ Wishart(n00), with n0 and Λ0 known;

where i = 1 for ψ ~ MCAR(a,Λ) and i = 1, 2 for ψ ~ MCAR(a1, a2,Λ)).

To avoid range restrictions in parameters ai, considering the transformations ρi = log(ai/(1−ai)) ∈ ℝ, then, the joint prior density for the parameters ϑ = (b, β, κ, ζ, ψ,Λ, ρ) can be written as

$π(ϑ)=π(b)π(β)π(ζ)π(ψ∣Λ,ρ)π(Λ)π(ρ)×exp(∑i=1Qζi)∏i=12exp(-ρi)(1+exp(-ρi))2.$

The joint posterior density is given by

$π(ϑ∣D)∝L(ϕ-1)exp{-12[σb-2∑i=0p1bi2+σβ-2∑i=1p2βi2+∑i=1Qexp(2ζi)σα2]+ψ′[Λ⊗(DW-aW)] ψ+log∣Λ⊗aW∣+n0-42log∣Λ∣-12tr(Λ0-1Λ)+∑i=1Qζi}π(ρ),$

where ϕ−1 = (b, β, ζ−1, U, V) and π(ρi) = 1 if ai ~ Uniform(0, 1) and π(ρi) = {1/B(18, 2)}{exp(17ρi)/(1 + exp(ρi))18} if ai ~ Beta(18, 2), where B(18, 2) = 17!/18! = 1/18.

This joint posterior density is analytically intractable. Thus, we based our inference on MCMC simulation methods. We can observed that the full conditional distributions for the parameters b, β, ζ, ψ and ρ do not present closed form; therefore, we will use the Metropolis-Hastings algorithm to generate posterior samples for these parameters. However, the Gibbs sampler algorithm is used to generate a posterior sample for the parameter Λ, because its full conditional distribution has a closed form. The full conditional distribution π(Λ|ϑ(−Λ), Dobs) is proportional to

$π(ψ∣Λ,a)π(Λ)∝ ∣Λ⊗DW-aW∣12exp(-12ψ′(DW-aW)ψ)∣Λ∣n0-42exp(-12tr(Λ0-1Λ))∝ ∣Λ∣I+n0-42exp{-12tr((Λ0-1+B)Λ)},$

where

$B=[tr(R1U (R1U)′)tr (R1U(R2V)′)tr(R2V (R1U)′)tr (R2V(R2V)′)].$

Thus, the full conditional distribution for Λ can be taken by the Wishart distribution with scale matrix $(Λ0-1+B)-1$ and degrees of freedom I + n0.

### 3.3. Model comparison criteria

There are several Bayesian criteria to compare competing models for a given data set and to select the one that best fits the data. One of the most used in applied works is based on the posterior mean of deviance and called the DIC. For a model, the statistic DIC is defined as

$DIC=d¯+pd,$

where = E[D(ϕ)], pd = E[D(ϕ)] − D[E(ϕ)] and D(ϕ) is the deviance function of the model defined by −2 log L(ϕ). L is the likelihood function of the model. Spiegelhalter et al. (2002) provide evidences that pd is a suitable measure of model complexity even in hierarchical settings, and thus, DIC is considered as a sensible generalization of the expected Akaike information criterion to hierarchical settings. The model, with the smallest value of DIC, is commonly taken as the preferred model describe the data set given.

### 3.4. Bayesian case influence diagnostics

Performing a sensitivity analysis is advisable since regression models are sensitive to underlying model assumptions. One of the most used ways of evaluating the influence of an observation in the fitted model is a case-deletion (Cook andWeisberg, 1982), where the effects are studied by completely removing cases from the analysis. This reasoning forms the basis of the Bayesian global influence methodology and makes it possible to determine which subjects might influence the analysis. The Bayesian case-deletion influence diagnostic measures for the joint posterior distribution based on the ψ-divergence (Peng and Dey, 1995; Weiss, 1996) is now introduced as follows.

Let Dψ(P, P(−i)) denote the ψ-divergence between P and P(−i), where P denotes the posterior distribution of ϑ for full data, and P(−i) denotes the posterior distribution of ϑ without the ith case. Therefore,

$Dψ(P,P(-i))=∫ψ(π(ϑ∣D(-i))π(ϑ∣D))π(ϑ∣D) dϑ,$

where ψ is a convex function with ψ(1) = 0. Several choices concerning the ψ are given by Dey and Birmiwal (1994). For example, ψ(z) = − log(z) defines the Kullback-Leibler (K-L) divergence, ψ(z) = (z − 1) log(z) gives J-distance (or the symmetric version of K-L divergence), ψ(z) = 0.5|z − 1| defines the variational distance (or L1 norm) and ψ(z) = (z − 1)2 defines the χ2-square divergence.

Let ϑ(1), …, ϑ(Q) be a size Q sample of , Dψ(P, P(−i)) can be calculated numerically by

$Dψ^(P,P(-i))=1Q∑q=1Qψ(CPO^iL(yi∣ϑ(q))),$

where $CPO^i=[(1/Q)∑q=1Q1/{L(yi∣ϑ(q))}]-1$ is the numerical approximation of the conditional predictive ordinate statistic of ith observation (Ibrahim et al., 2001).

Note that Dψ(P, P(−i)) can be interpreted as the ψ-divergence of the effect of deleting the i-th case from the full data on the joint posterior distribution of ϑ. As pointed by Peng and Dey (1995), Weiss (1996), and Cancho et al. (2010), it may be difficult for a practitioner to judge the cutoff point of the divergence measure to determine if a small subset of observations is influential. In this context, we will use the proposal given by Peng and Dey (1995) and Weiss (1996) by considering a biased coin, which has success probability p. Then the ψ-divergence between the biased and an unbiased coin is

$Dψ(f0,f1)=∫ψ(f0(x)f1(x))f1(x)dx,$

where f0(x) = px(1− p)1−x and f1(x) = 0.5, x = 0, 1. Now if Dψ( f0, f1) = dψ(p), then it would be easy to check if dψ satisfies the following equation

$dψ(p)=ψ(2p)+ψ(2(1-p))2.$

It is not difficult to see for the divergence measures considered that dψ increases as p moves away from 0.5. In addition, dψ(p) is symmetric about p = 0.5 and dψ, achieves its minimum at p = 0.5. At this point, dψ(0.5) = 0, and f0 = f1. Therefore, if we consider p > 0.90 (or p ≤ 0.10) as a strong bias in a coin, then dK-L(0.90) = 0.51, dJ(0.90) = 0.88, dL1(0.90) = 0.4, and dχ2(0.90) = 0.64. This equation implies that ith case is influential when dL1 > 0.4 or dχ2 > 0.64. Thus, if we use the K-L divergence, we can consider an influential observation when dK-L > 0.51. Similarly, using the J-distance, an observation which implies dJ > 0.88 can be considered influential.

4. Simulation study

This section presents simulation studies for a cure rate model under FA with the dependent assumption of examining their performances. The interval-censored survival times (ti jL, ti jR, δi j) with the cure fraction under the FA are generated similarly to what was introduced by Yau and Ng (2001) with some modifications.

We generate the latent Geometric variable Mi j, which denotes the initial number of competing causes regarding the event, with parameter p0i j = [1 + exp(−(b0 + b1xi j) + vi)]−1 for the jth individual in the ith region, j = 1, …, ni, i = 1, …, I, where the covariate xi j is generated from the Bernoulli distribution with parameter 0.5. Interval-censored data (ti jL, ti jR, δi j) are then generated as:

• If Mi j = 0, then let ti j = ti jL from the Exponential distribution with hazard rate 10 and let the censoring indicator δi j = 0.

• If Mi j > 0, then we generate Mi j latent Weibull variables with parameter α = 0.30 and λi j = (βxi j + ui) with β = −0.15. Let ti j takes the lowest generated variable in case of generating the variables of the model under FA. The censoring times ci j were sampled from Uniform U(0, τ), where τ > 0 was set in order to control the proportion of censored observations to be approximately 40% on average and, the indicator variables will be of δi j = 1 if ti jci j and δi j = 0, otherwise.

• For δi j = 0, let 0 < ti jL < ti jR = ∞.

• For δi j = 1, we create leni j from the distribution U(0.2, 0.7) and li j from U(0, 0.01). Then, from (0, li j], (li j, li j + leni j], …, (li j + kleni j,∞], k = 1, 2, …, the interval (ti jL,ti jR] which satisfies ti jL < ti jti jR is chosen.

In this study, we consider I = 5 regions (Zip) with the corresponding adjacent matrix

$[0010000001100100010001000].$

Random effects Ui and Ui are generated from the Normal distribution with mean 0 and precision matrix Λ ⊗ (DWaW), where W is the standardized adjacent matrix so that each of its rows sum one, DW = Diag(1, 1, 2, 1, 1) is a diagonal matrix, and then we fixed a = 0.9 and Λ = Diag(4, 4), i.e., we fixed Λ11 = 4, Λ22 = 4, and Λ12 = Λ21 = 0. By considering 100 individuals in the simulation studies. The corresponding zip codes for each individual were distributed using sample with replace, thus the number of individuals in each region ni, i = 1, …, 5 are varied; therefore, these five regions could present different numbers of individuals with $∑i=15ni=100$. Thus, we have sample size n = 100 and we fixed the parameters b0 = −1.50, b1 = −0.50, β = −0.15, α = 0.30. In simulations, we used to consider around 40% of the censored data for each generated sample with 500 repeated samples simulated for each model. The priors for the parameters b0, b1, β1, and α used in the studies are b0 ~ N(0, 32), b1 ~ N(0, 32), β1 ~ N(0, 32), and α ~ N(0, 102)I(0,∞).

For each generated sample, we simulate one chain of size 10,000 for each parameter, disregarding the first 1,000 iterations to eliminate the effect of the initial values and avoid correlation problems and thinning to every third iteration, thus obtaining an effective sample of size 3,000 upon which the posterior is based on. To evaluate the performance of the parameter estimates, the average bias (Bias), standard deviation (SD) of the estimate, average SD (SDs mean), and the mean square error (MSE) are calculated for a cure rate model under the FA (Table 1).

We can note that the bias and the MSE of the parameter Λ12 are larger than other parameters. The estimator of the Λ12 presents a negative biases; however, the biases and MSEs are always near zero. Moreover, the simulation results for the cure model considering the prior 1 are close to those obtained considering the prior 2.

### Influence of outlying observations

A goal of this study is to show the need for robust models to deal with the presence of outliers in data. We consider two cases for perturbation with the parameter’s values and the setup the same as in the simulation studies; therefore, four data sets of size 100 were generated from the cure model under the FA with dependent spatial frailties.

We selected the cases 18 and 80 for perturbation. To create influential observations concerning the data set, we choose one or two of these selected cases and perturbed the response variable as: $tkL˜=tkL+10SDL$ and $tkR˜=tkR+10SDL$, for k = 18 and 80, where SDL is the SD of the ti jL. Note that this perturbation method will not change the time interval of the perturbation candidate. Here, we consider four setups in the study. Setup A: original dataset, without outliers, Setup B: data with outlier 18, Setup C: data with outlier 80, and Setup D: data with outliers 18 and 80. The MCMC computations were similar to those in the last section; in addition, Gewekes convergence diagnostic proposed by Geweke (1992) was used to monitor the convergence of the Gibbs samples.

Table 2 reports the posterior mean, the SD, the bias and the MSE of the parameters of WCRM-FA. We can note that the estimative of parameter Λ11 creasing in the perturbation cases when prior 1 is used. However, considering prior 2 for the parameters, the estimative of all parameters of cases B, C, and D are very closed the case A, which means that the parameters are not sensitive to perturbations.

For each simulated data set, the four divergence measures (dKL, dJ, dL1,dχ2) of the perturbed cases and the DIC values for theWeibull and the proportional hazards (PH) cure rate models were calculated (Table 3).

We note that all measures provide larger ψ-divergence measures when compared to the non-perturbed setup (setup A) and that the difference between the measures of the perturbed cases and non-perturbed cases is more clear for the PH cure rate models than the WCRMs. Furthermore, we observe that the values of the measures obtained from the cure models, whether the parameters have the prior 1 or 2, are similar.

To better show present the results, we plot the J-distance measure from the fitted models considering the prior 1 for the parameters. The Figure 1 presents the divergence measures before the perturbation (setup A), the model indicates the absence of outline observations, and after perturbation observations (setups B, C, and D). Note that outline observation 18 cannot be easy detected for the cure rate model under the FA.

5. Application

Recall the interval-censored smoking cessation data briefly presented in the introduction section. In the smoking cessation study, all patients (smokers) were randomized into either a SI group, or a UC group that received no special anti-SI. The treatments of the SI program were realized in Rochester city localizing in the center of the maps. The program detail can be found in Murray et al. (1998). Here, each patient was observed annually over five years of monitoring. Our event of interest concerns if they will relapse (resume smoking). If a smoker starts smoking again after an initial attempt to quit, then only an approximate one-year time interval was observed from the previous observation to the current observation. Thus, the relapse times are interval-censored. In this analysis, we limit our attention to those patients who are known to have quit smoking at least once during the study period and who have an identifiable Minnesota zip code of residence. Thus, the data consists of a total of 223 patients who reside in 51 zip codes in the Southeastern corner of Minnesota, among them there are 65 patients that have undergone relapse, which implies the empirical cure rate is approximately 71%. Moreover, Ma and Xiang (2013) have also confirmed the existence of non-eligible cure fraction in the population.

We fitted cure rate model under first and LA, considering the different spatial frailties in the models to the data set. The prior distributions for the parameters b, β, and α are:

• bj ~ N(0, 100), j = 0, …, 4;

• βj ~ N(0, 100), j = 1, …, 4;

• α ~ N(0, 100)I(0,∞).

Because of the high computational cost, we implement the MCMC algorithms in the C programming language; subsequently, the results were analyzed in the R language (R Development Core Team, 2010) through the “coda” package (Plummer et al., 2005). All of our MCMC algorithms ran a total of 60,000 iterations discarding the first 20,000 realizations as burn-in and thinning to every fifth iteration. Posterior results are then based on 8,000 realizations of the Markov chain. Our Metropolis acceptance rate for these parameters ranged from 25% to 50%. The convergence was checked using the Geweke diagnostic which did not indicate a lack of convergence. The models are compared using DIC criterion.

Table 4 provides the DIC scores for a variety of effects of the Weibull cure model under FA and LA. The DIC scores of the Model 1 and 5 show the best models despite DIC values that are close to each other. We also can note that the cure rate models under the FA are more adequate than models under the LA.

For the comparison with the models proposed by Carlin and Banerjee (2003), we consider the same prior distributions for the parameters b and β as considered by these authors. Table 5 reports the DIC scores for a variety of effects for the WCRM. We observe that the DIC scores in Table 5 are smaller than the values in Table 4. However, the DIC scores for the cure rate model under FA are very close each other and indicates that these models are equivalent. Moreover, similarly to the previous case, the DIC values of the cure rate model under FA are smaller than the cure rate model under LA. Comparing the obtained DIC scores in Table 4 and 5 with the DIC scores presented in the study conducted by Carlin and Banerjee (2003), where they proposed the mixture cure model with the spatial frailty, we can conclude that all our models are more adequate since the DIC scores are smaller. Here, we select the Model 15 as our working model.

Table 6 presents the posterior summary of the parameters of the Model 15. We note that only the parameters b0, b2, b3, and β3 are significant. In the cure rate, the negative value of b3 means that the individuals with higher levels of cigarette consumption have a lower probability to quit smoking, while the positive value of b2 implies the individuals with special intervention have higher probability to quit smoking than those with UC.

The survival model shows that the special intervention and the number of cigarettes smoked per day have negative effects on the hazard rate of the relapse time; therefore, individuals with special intervention do not present lower hazard rates for the relapse time when compared to those who attend UC. However, the individuals with a higher level of cigarette consumption do not present high hazard rates.

The estimated $SD Σ111/2$ of random spatial effects in the survival model is 0.4113, and the estimated $SD Σ221/2$ of the random spatial effects in the cure rate is 0.4268, which indicates that there is a considerable heterogeneity among clusters. It is also observed that there are no correlations between the spatial effects U and V.

The Figure 2 maps the posterior means of the frailties U and V in the Model 15. For the frailties U of which the high value presents the high relapse rate, we note that the city of Owatonna and some southern cities have higher values with individuals in these regions indicating higher relapse rates than others. However, the city of Rochester (in the central region) suggests slightly better than avenge cessation behavior, which can also be observed by the frailties V. Note that the high value of V presents the high cure probability.

The Figure 3 maps the posterior SD of the frailties U and V corresponding to the posterior means mapped in Figure 2. We note the posterior SD of the frailties U and V have approximated values. Both maps show that the central region cities have lower values and that some periphery cites have higher values.

Figure 4 presents the estimates of ψ-divergence measures, which were obtained by the posterior sample of the parameters of the model that are used to detect possible influential observations in the posterior distribution of the parameters of Model 15. It shows some individuals can be influential observations that can be detected by divergence measures. Here, we will analyze individuals 72, 138, 151, and 199 were detected by the J-distance and the χ2-divergence. Table 7 presents information on them, so that we can note that these four individuals had special interventions and did not consume high amounts of cigarettes per day, while the individuals 72, 138, and 151 had relapse, but the individual 199 did not. To reveal the impact of these possible influential observations on the parameter estimates and inferences, we removed such observations, refitting the models. We also calculated the relative variations (RV) for the posterior mean of the parameters. The RV is defined by RV = (ϑ̂d,−Iϑ̂d)/ϑ̂d, for all d, where I denotes a set of influential observations, d is the index of the parameters, ϑ̂d,−I denotes the posterior mean of ϑd,−I, after the set of observations I was removed. In this case, we have I = {72, 138, 151, 199}.

The posterior summaries of the parameters for the readjust Model 15 and RV for the posterior mean of the parameters are presented in Table 8. We can note that only the values of RV for the posterior means of the parameters Λ12 and ∑12 are more than one, but they still have posterior means near zero, and others parameters have posterior means near the obtained values for the completed data set. In this case, there are no inferential changes after removing the observations.

6. Conclusions

This paper described an approach to extend the proportional odds cure models to allow to spatial correlations by including spatial frailties in the interval-censored data setting. We used MCMC methods through Bayesian inference and the DIC for the model comparison. The results of the application show that the parametric cure model with frailties under the FA scheme have better fittings. A comparison of the proposed models with the model introduced by Carlin and Banerjee (2003) indicated that our model is more adequate. Moreover, the both models proposed are not sensitive to influence observations, which can be observed by the influence diagnostic analysis in the application. The interpretation of the covariates was easy due to the parametrization of the models considered in terms of the cure rate. The MCAR prior can also be used even if frailties do not present or do present low correlations. We assume Mi j as Geometric distributed; however, we believe that the methodology presented here may be theoretically extended by considering other distributions, such as Negative Binomial and Power Series.

Acknowledgements

The research is partially supported by the Brazilian organizations CAPES (via Science without Borders-PDSE), CNPq and FAPESP.

Figures
Fig. 1. Index plots of J-distance measure from the fitted of the Weibull cure rate model-first activation. From left to right, the data sets of the figures correspond to the setup A, B, C, and D.
Fig. 2. Maps of posterior means for frailties (left panel) and (right panel) in Model 15.
Fig. 3. Maps of posterior standard deviation for frailties (left panel) and (right panel) in Model 15.
Fig. 4. Estimates of ψ-divergence measures for Model 15. K-L = Kullback-Leibler.
TABLES

### Table 1

Simulation results for the Weibull cure model under the first activation with depended spatial frailties

ParameterTrue valueEstimate meanSD of the estimateBiasMSESDs mean
Prior 1: ψ ~ MCAR(a,Λ), a ~ Beta(18, 2), Λ0 ~ Wishart(2, Diag(0.9, 1))

b0−1.50−1.48120.07080.01880.00540.2685
b1−0.50−0.52090.1280−0.02090.01680.2485
β−0.15−0.13600.04660.01400.00240.1915
α0.300.19300.0521−0.10700.01420.0677
Λ114.004.00620.15970.00620.02552.4728
Λ224.004.01200.19180.01200.03692.6151
Λ120.00−0.45410.1343−0.45410.22421.9196
a0.900.90010.00160.00010.00000.0653

Prior 2: ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Beta(18, 2), Λ0 ~ Wishart(2, Diag(0.9, 1))

b0−1.50−1.49020.07010.00980.00500.2583
b1−0.50−0.53760.1330−0.03760.01910.2227
β−0.15−0.12950.04930.02050.00280.1870
α0.300.18630.0443−0.11370.01490.0536
Λ114.004.16380.16760.16380.05492.5070
Λ224.004.26570.18190.26570.10362.6919
Λ120.00−0.58090.1472−0.58090.35911.9647
a10.900.89990.0015−0.00020.00000.0655
a20.900.90020.00150.00020.00000.0654

SD = standard deviation; Bias = average bias; MSE = mean square error; SDs mean = average SD; MCAR = multivariate conditionally autoregressive.

### Table 2

Simulation results of the perturbed cases for the Weibull cure rate models under the first activation

SetupPerturbed casePrior 1Prior 2

ParametersMeanSDBiasMSEParametersMeanSDBiasMSE
ANoneb0−1.4820.2780.0180.000b0−1.3410.2770.1590.025
b1−0.8590.259−0.3590.129b1−0.5540.248−0.0540.003
β−0.0360.1930.1140.013β−0.1050.1970.0450.002
α0.2280.068−0.0720.005α0.1120.072−0.1880.035
Λ114.1722.5330.1720.030Λ114.0522.4860.0520.003
Λ224.1522.6930.1520.023Λ224.0122.6360.0120.000
Λ12−0.5021.945−0.5020.252Λ12−0.5401.905−0.5400.291
a0.9020.0640.0020.000a10.8990.068−0.0010.000
a20.8990.067−0.0010.000

B{18}b0−1.5260.263−0.0260.001b0−1.5720.275−0.0720.005
b1−0.4390.2530.0610.004b1−0.6220.260−0.1220.015
β−0.1550.189−0.0050.000β−0.0990.1900.0510.003
α0.2310.061−0.0690.005α0.3320.0730.0320.001
Λ114.0572.4330.0570.003Λ114.0562.4940.0560.003
Λ223.9472.630−0.0530.003Λ223.9132.561−0.0870.008
Λ12−0.4971.904−0.4970.247Λ12−0.4061.991−0.4060.164
a0.9000.0650.0000.000a10.8990.065−0.0010.000
a20.9000.0650.0000.000

C{80}b0−1.5100.275−0.0100.000b0−1.5240.263−0.0240.001
b1−0.5860.257−0.0860.007b1−0.6790.264−0.1790.032
β−0.1390.1940.0110.000β−0.0830.1930.0670.005
α0.2550.082−0.0450.002α0.2350.078−0.0650.004
Λ113.8322.406−0.1680.028Λ114.0182.3860.0180.000
Λ223.5352.423−0.4650.216Λ223.9192.649−0.0810.007
Λ12−0.1151.835−0.1150.013Λ12−0.3021.969−0.3020.091
a0.9000.0630.0000.000a10.8990.064−0.0010.000
a20.9000.0680.0000.000

D{18, 80}b0−1.4600.2590.0400.002b0−1.5990.266−0.0990.010
b1−0.3480.2480.1520.023b1−0.3160.2460.1840.034
β−0.2030.189−0.0530.003β−0.2080.191−0.0580.003
α0.1870.059−0.1130.013α0.2530.061−0.0470.002
Λ114.2052.5690.2050.042Λ114.2592.5410.2590.067
Λ224.0862.6340.0860.007Λ224.0022.5490.0020.000
Λ12−0.5151.971−0.5150.265Λ12−0.5881.946−0.5880.346
a0.9000.0640.0000.000a20.9010.0650.0010.000
a20.8980.068−0.0020.000

SD = standard deviation; Bias = average bias; MSE = mean square error.

### Table 3

ψ-divergence measures of the perturbed cases and the values of DIC for simulated data sets

ActivationPriorSetupPerctausrebedNo. of casedKLdJdL1dχ2DIC
First1ANone180.0060.0110.0420.012142.754
800.0300.0600.0970.067
B{18}180.0800.1710.1590.238164.476
C{80}800.2460.6020.2772.128153.082
D{18,80}180.0690.1430.1470.181185.709
800.2820.6770.3041.742

2ANone180.0070.0140.0460.014140.186
800.0330.0670.1020.075
B{18}180.0360.0750.1060.084164.446
C{80}800.2940.9400.28814.544149.760
D{18,80}180.0620.1310.1380.176184.934
800.1200.2690.1900.498

DIC = deviance information criterion.

### Table 4

Real dataset: Bayesian criteria for the fitted models

Cure rate model under first activationCriteria

ModelPriorsDICpd
1U ~ CAR(ρ1), V ~ CAR(ρ2), ρ1, ρ2 ~ InvGamma(0.01, 0.01)416.391711.5105
2ψ ~ MCAR(a,Λ), a ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))416.969111.6971
3ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))417.494212.1639
4ψ ~ MCAR(a,Λ), a ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))417.156012.6639
5ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))416.670011.7184

Cure rate model under last activationCriteria

ModelPriorsDICpd

6U ~ CAR(ρ1), V ~ CAR(ρ2), ρ1, ρ2 ~ InvGamma(0.01, 0.01)419.272311.8914
7ψ ~ MCAR(a,Λ), a ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))419.158512.2811
8ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))419.450912.6852
9ψ ~ MCAR(a,Λ), a ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))417.891713.6477
10ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))418.643213.9228

DIC = deviance information criterion.

### Table 5

Real dataset: Bayesian criteria for the fitted models

Cure rate model under first activationCriteria

ModelPriorsDICpd
11U ~ CAR(ρ1), V ~ CAR(ρ2), ρ1, ρ2 ~ InvGamma(0.01, 0.01)414.37518.2277
12ψ ~ MCAR(a,Λ), a ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))414.762010.9257
13ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1)) 1414.748410.8353
14ψ ~ MCAR(a,Λ), a ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))414.348310.9108
15ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))414.492510.9244

Cure rate model under last activationCriteria

ModelPriorsDICpd

16U ~ CAR(ρ1), V ~ CAR(ρ2), ρ1, ρ2 ~ InvGamma(0.01, 0.01)418.13369.4144
17ψ ~ MCAR(a,Λ), a ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))417.271711.7087
18ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Uniform(0, 1), Λ ~ Wishart(2, Diag(0.1, 0.1))416.862011.4699
19ψ ~ MCAR(a,Λ), a ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))416.897911.6343
20ψ ~ MCAR(a1, a2,Λ), a1, a2 ~ Beta(18, 2), Λ ~ Wishart(2, Diag(0.1, 0.1))416.878211.5676

DIC = deviance information criterion.

### Table 6

Posterior summaries of the parameter of the Model 15 for the smoking cessation data

ParameterSurvival modelCure rate

MeanSD2.5%97.5%MeanSD2.5%97.5%
Interceptb01.37360.57890.26002.5442
Sex (male = 0)β1−0.15620.4551−1.06260.6992b1−0.50960.3718−1.27110.2052
SI/UC (UC = 0)β20.84270.5191−0.16591.8703b20.86010.43100.06711.7192
Cigarettes per dayβ3−0.11480.0378−0.1809−0.0322b3−0.07280.0290−0.1345−0.0201
Duration as smokerβ4−0.02460.0343−0.10030.0345b40.01970.0230−0.02640.0651
α2.40970.30731.81133.0065
a1(au)0.89680.06760.72610.9874
a2(av)0.89940.06700.73700.9881

Λ112.67690.63771.55434.0673
Λ222.57430.64131.49243.9718
Λ12−0.01040.4625−0.93210.8919

110.41130.10750.25310.6754
220.42980.12360.25610.7298
12/(∑1122)1/20.00350.1828−0.36550.3559

Where Λi j is the element of precision matrix Λ in position (i, j), and ∑i j is the element of matrix Λ = Λ−1 in position (i, j), this ∑11 is the spatial variance component of U and ∑22 is the spatial variance component of V, ∑12/(∑1122)1/2 denote their correlation. SD = standard deviation; SI = smoking intervention group; UC = usual care group.

### Table 7

Possible influential observations are detected by four divergence measures

ObservationsSexDurationInterventionNo. of cigarettesRelapseTime intervalZip
720201251(3.159, 3.929)55987
1381251201(2.998, 3.992)55021
1510391101(0.923, 3.962)55057
1990220201(3.885, 5.013)55904

### Table 8

Posterior summaries of the parameter of the Model 15 and relative variations adjusted for the smoking cessation data without detected individuals 72, 138, 151, and 199

ParameterSurvival modelCure rate

MeanSD2.5%97.5%MeanSD2.5%97.5%
Interceptb0 (−0.0028)1.36970.59810.20412.5265

Sex (male = 0)β1−0.2101 (−0.3455)0.4326−1.08480.6344b1−0.5657 (0.1101)0.3425−1.25110.0902

SI/UC (UC = 0)β21.0138 (0.2031)0.50710.04862.0473b2 (0.0491)0.90230.40020.14061.7522

Cigarettes per dayβ3−0.1048 (−0.0873)0.0377−0.1800−0.0304b3 (−0.1602)−0.06110.0257−0.1221−0.0183

Duration as smokerβ4−0.0308 (0.2518)0.0361−0.10730.0375b40.0184 (−0.0656)0.0230−0.02800.0629

α2.6474 (0.0987)0.35151.97923.3582

a1 (au)0.9006 (0.0043)0.06650.73880.9880

a2 (av)0.9002 (0.0009)0.06500.73740.9864

Λ112.6749 (−0.0008)0.64971.57694.0686

Λ222.5717 (−0.0010)0.63641.51783.9480

Λ120.0113 (−2.0812)0.4709−0.92570.9130

110.4122 (0.0023)0.10860.25150.6701

220.4307 (0.0021)0.12080.25750.7207

12−0.0054 (−2.5601)0.1839−0.36370.3568

SD = standard deviation; SI = smoking intervention group; UC = usual care group.

References
1. Banerjee, S, and Carlin, BP (2004). Parametric spatial cure rate models for interval-censored time-to-relapse data. Biometrics. 60, 268-275.
2. Bennett, S (1983). Analysis of survival data by the proportional odds model. Statistics in Medicine. 2, 273-277.
3. Berkson, J, and Gage, RP (1952). Survival curve for cancer patients following treatment. Journal of the American Statistical Association. 47, 501-515.
4. Besag, J (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society, Series B (Methodological). 36, 192-236.
5. Cancho, VG, Ortega, EMM, and Paula, GA (2010). On estimation and influence diagnostics for log-Birnbaum-Saunders Student-t regression models: full Bayesian analysis. Journal of Statistical Planning and Inference. 140, 2486-2496.
6. Carlin, BP, and Banerjee, S (2003). Hierarchical multivariate CAR models for spatio-temporally correlated survival data. Bayesian Statistics. 7, 45-63.
7. Chen, MH, Ibrahim, JG, and Sinha, D (1999). A new Bayesian model for survival data with a surviving fraction. Journal of the American Statistical Association. 94, 909-919.
8. Cook, RD, and Weisberg, S (1982). Residuals and Influence in Regression. FL: Chapman & Hall/CRC
9. Collett, D (1994). Modelling survival data. Modelling Survival Data in Medical Research, Collett, D, ed. New York: Springer, pp. 53-60
10. Cooner, F, Banerjee, S, Carlin, BP, and Sinha, D (2007). Flexible cure rate modeling under latent activation schemes. Journal of the American Statistical Association. 102, 560-572.
11. Dey, DK, and Birmiwal, LR (1994). Robust Bayesian analysis using divergence measures. Statistics & Probability Letters. 20, 287-294.
12. Finkelstein, DM (1986). A proportional Hazards model for interval-censored failure time data. Biometrics. 42, 845-854.
13. Gamerman, D, and Lopes, HF (2006). Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. Boca Raton: Chapman & Hall/CRC
14. Gelfand, AE, and Vounatsou, P (2003). Proper multivariate conditional autoregressive models for spatial data analysis. Biostatistics. 4, 11-15.
15. Geweke, J (1992). Evaluating the accuracy of sampling-based approaches to calculating posterior moments. Bayesian Statistics 4, Bernardo, JM, Berger, J, Dawid, AP, and Smith, JFM, ed. Oxford: Clarendon Press, pp. 169-193
16. Gu, Y, Sinha, D, and Banerjee, S (2011). Analysis of cure rate survival data under proportional odds model. Lifetime Data Analysis. 17, 123-134.
17. Hu, T, and Xiang, L (2013). Efficient estimation for semiparametric cure models with interval-censored data. Journal of Multivariate Analysis. 121, 139-151.
18. Ibrahim, JG, Chen, MH, and Sinha, D (2001). Bayesian Survival Analysis. New York: Springer
19. Lindsey, JC, and Ryan, LM (1998). Methods for interval-censored data. Statistics in Medicine. 17, 219-238.
20. Ma, X, and Xiang, L (2013). Testing for the presence of a cure fraction in clustered interval-censored survival data. Australian & New Zealand Journal of Statistics. 55, 173-190.
21. Murray, RP, Anthonisen, NR, and Connett, JE (1998). Effects of multiple attempts to quit smoking and relapses to smoking on pulmonary function. Journal of Clinical Epidemiology. 51, 1317-1326.
22. Pan, C, Cai, B, Wang, L, and Lin, X (2014). Bayesian semiparametric model for spatially correlated interval-censored survival data. Computational Statistics & Data Analysis. 74, 198-208.
23. Peng, F, and Dey, DK (1995). Bayesian analysis of outlier problems using divergence measures. Canadian Journal of Statistics. 23, 199-213.
24. Peto, R (1973). Experimental survival curves for interval-censored data. Journal of the Royal Statistical Society Series C (Applied Statistics). 22, 86-91.
25. Plummer, M, Best, N, Cowles, K, and Vines, K (2005). Output analysis and diagnostics for MCMC: R package version 0.10–3.Retrieved November 2, 2017, from: http://cran.rproject.org
26. , (2010). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing
27. Rodrigues, J, Cancho, VG, de Castro, M, and Louzada-Neto, F (2009). On the unification of long-term survival models. Statistics & Probability Letters. 79, 753-759.
28. Rücker, G, and Messerer, D (1988). Remission duration: an example of interval-censored observations. Statistics in Medicine. 7, 1139-1145.
29. Spiegelhalter, DJ, Best, NG, Carlin, BP, and van der Linde, A (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B (Statistical Methodology). 64, 583-639.
30. Sun, X, and Chen, C (2010). Comparison of Finkelstein’s method with the conventional approach for interval-censored data analysis. Statistics in Biopharmaceutical Research. 2, 97-108.
31. Weiss, R (1996). An approach to Bayesian sensitivity analysis. Journal of the Royal Statistical Society Series B (Methodological). 58, 739-750.
32. Xiang, L, Ma, X, and Yau, KK (2011). Mixture cure model with random effects for clustered interval-censored survival data. Statistics in Medicine. 30, 995-1006.
33. Yakovlev, AY, and Tsodikov, AD (1996). Stochastic Models of Tumor Latency and Their Biostatistical Applications. Singapore: World Scientific
34. Yau, KKW, and Ng, ASK (2001). Long-term survivor mixture model with random effects: application to a multi-centre clinical trial of carcinoma. Statistics in Medicine. 20, 1591-1607.
35. Zeng, D, Yin, G, and Ibrahim, JG (2006). Semiparametric transformation models for survival data with a cure fraction. Journal of the American Statistical Association. 101, 670-684.