spBayesSurv: Fitting Bayesian Spatial Survival Models Using R
JJSS
Journal of Statistical Software
MMMMMM YYYY, Volume VV, Issue II. doi: 10.18637/jss.v000.i00 spBayesSurv : Fitting Bayesian Spatial SurvivalModels Using R Haiming Zhou
Northern Illinois University
Timothy Hanson
Medtronic Inc.
Jiajia Zhang
University of South Carolina
Abstract
Spatial survival analysis has received a great deal of attention over the last 20 yearsdue to the important role that geographical information can play in predicting survival.This paper provides an introduction to a set of programs for implementing some Bayesianspatial survival models in R using the package spBayesSurv . The function survregbayes includes the three most commonly-used semiparametric models: proportional hazards,proportional odds, and accelerated failure time. All manner of censored survival timesare simultaneously accommodated including uncensored, interval censored, current-status,left and right censored, and mixtures of these. Left-truncated data are also accommodated.Time-dependent covariates are allowed under the piecewise constant assumption. Bothgeoreferenced and areally observed spatial locations are handled via frailties. Model fitis assessed with conditional Cox-Snell residual plots, and model choice is carried out viathe log pseudo marginal likelihood, the deviance information criterion and the Watanabe-Akaike information criterion. The accelerated failure time frailty model with a covariate-dependent baseline is included in the function frailtyGAFT . In addition, the packagealso provides two marginal survival models: proportional hazards and linear dependentDirichlet process mixture, where the spatial dependence is modeled via spatial copulas.Note that the package can also handle non-spatial data using non-spatial versions ofaforementioned models. Keywords : Bayesian nonparametric, survival analysis, spatial dependence, semiparametricmodels, parametric models.
1. Introduction
Spatial location plays a key role in survival prediction, serving as a proxy for unmeasuredregional characteristics such as socioeconomic status, access to health care, pollution, etc. Lit-erature on the spatial analysis of survival data has flourished over the last decade, including a r X i v : . [ s t a t . C O ] A p r spBayesSurv, version 1.1.3 the study of leukemia survival (Henderson, Shimakura, and Gorst 2002), childhood mortal-ity (Kneib 2006), asthma (Li and Lin 2006), breast cancer (Banerjee and Dey 2005; Zhou,Hanson, Jara, and Zhang 2015a), political event processes (Darmofal 2009), prostate cancer(Wang, Zhang, and Lawson 2012; Zhou, Hanson, and Zhang 2017), pine trees (Li, Hong,Thapa, and Burkhart 2015), threatened frogs (Zhou, Hanson, and Knapp 2015b), health andpharmaceutical firms (Arbia, Espa, Giuliani, and Micciolo 2016), emergency service responsetimes (Taylor 2017), and many others.Here we introduce the spBayesSurv (Zhou and Hanson 2018) package for fitting various sur-vival models to spatially-referenced survival data. Note that all models included in thispackage can also be fit without spatial information, including nonparametric models as wellas semiparametric proportional hazards (PH), proportional odds (PO), and accelerated failuretime (AFT) models. The model parameters and statistical inference are carried out via self-tuning adaptive Markov chain Monte Carlo (MCMC) methods; no manual tuning is needed.The R syntax is essentially the same as for existing R survival (Therneau 2015) functions.Sensible, well-tested default priors are used throughout, however, the user can easily imple-ment informative priors if such information is available. The primary goal of this paper isto introduce spBayesSurv and provide extensive examples of its use. Comparisons to othermodels and R packages can be found in Zhou et al. (2015b), Zhou et al. (2017), and Zhou andHanson (2017).Section 2 discusses spBayesSurv ’s implementation of PH, PO, and AFT frailty models forgeoreferenced (e.g., latitude and longitude are recorded) and areally-referenced (e.g., countyof residence recorded) spatial survival data; the functions also work very well for exchangeableor no frailties. The models are centered at a parametric family through a novel transformedBernstein polynomial prior and the centering family can be tested versus the Bernstein ex-tension via Bayes factors. All manner of censoring is accommodated as well as left-truncateddata; left-truncation also allows for the inclusion of time-dependent covariates. The LPML,DIC and WAIC statistics are available for model selection; spike-and-slab variable selectionis also implemented.In Section 3, a generalized AFT model is implemented allowing for continuous stratification .That is, the baseline survival function is itself a function of covariates: baseline survivalchanges smoothly as a function of continuous predictors; for categorical predictors the usualstratified AFT model is obtained. Note that even for the usual stratified semiparametric AFTmodel with one discrete predictor (e.g., clinic) it is extremely difficult to obtain inference usingfrequentist approaches; see Chiou, Kang, and Yan (2015) for a recent development. Themodel fit in spBayesSurv actually extends discrete stratification to continuous covariates,allowing for very general models to be fit. The generalized AFT model includes the easycomputation of Bayes factors for determining which covariates affect baseline survival andwhether a parametric baseline is adequate.Finally, Section 4 offers a spatial implementation of the completely nonparametric lineardependent Dirichlet process mixture (LDDPM) model of De Iorio, Johnson, M¨uller, andRosner (2009) for georeferenced data. The LDDPM does not have one simple “linear predictor”as do the models in Sections 2 and 3, and therefore a marginal copula approach was takento incorporate spatial dependence. A piecewise-constant baseline hazard PH model is alsoimplemented via spatial copula for comparison purposes, i.e., a Bayesian version of the modelpresented in Li and Lin (2006). Section 5 concludes the paper with a discussion. ournal of Statistical Software R packages for implementing survival models, there are only ahandful of that allow the inclusion of spatial information and these focus almost exclusivelyon variants of the PH model. BayesX (Belitz, Brezger, Klein, Kneib, Lang, and Umlauf 2015)is an immensely powerful standalone program for fitting various generalized additive mixedmodels, including both georeferenced and areally-referenced frailties in the PH model. Thepackage
R2BayesX (Umlauf, Adler, Kneib, Lang, and Zeileis 2015) interfaces
BayesX with R ,but does not appear to include the full functionality of BayesX , e.g., a Bayesian approach forinterval-censored data is not included.
BayesX uses Gaussian Markov random fields for discretespatial data. For georeferenced frailties
BayesX uses what have been termed “Matern splines,”first introduced in an applied context by Kammann and Wand (2003). Several authors haveused this approach including Kneib (2006), Hennerfeind, Brezger, and Fahrmeir (2006), andKneib and Fahrmeir (2007). This approximation was termed a “predictive process” and givena more formal treatment by Banerjee, Gelfand, Finley, and Sang (2008) and Finley, Sang,Banerjee, and Gelfand (2009). The spBayesSurv package utilizes the full-scale approximation(FSA) of Sang and Huang (2012) which extends the predictive process to capture both thelarge and small spatial scales; see Section 2.1.4.The package spatsurv (Taylor and Rowlingson 2017) includes an implementation of PH allow-ing for georeferenced Gaussian process frailties. The frailty process is approximated on a finegrid and the covariance matrix inverted via the discrete Fourier transform on block circulantmatrices; see Taylor (2015) for details. Taylor’s approach vastly improves computation timeover a fully-specified Gaussian process. The package mgcv (Wood 2017) also fits a spatialPH model by including a spatial term through various smoothers such as thin plate spline,Duchon spline and Gaussian process. All the three aforementioned R packages focus on thePH model, whereas the spBayesSurv includes several other spatial frailty models and twomarginal copula models (Zhou et al. m distinct spatial locations s , . . . , s m . Let t ij be a random event time associated with the j th subject in s i and x ij be a related p -dimensional vector of covariates, i = 1 , . . . , m, j = 1 , . . . , n i . Then n = (cid:80) mi =1 n i is the totalnumber of subjects under consideration. Assume the survival time t ij lies in the interval( a ij , b ij ), 0 ≤ a ij ≤ b ij ≤ ∞ . Here left censored data are of the form (0 , b ij ), right censored( a ij , ∞ ), interval censored ( a ij , b ij ) and uncensored values simply have a ij = b ij , i.e., we define( x, x ) = { x } . Therefore, the observed data will be D = { ( a ij , b ij , x ij , s i ); i = 1 , . . . , m, j =1 , . . . , n i } . For areally-observed outcomes, e.g., county-level, there is typically replication(i.e., n i > R version 3.3.3 under the platform x86 64-apple-darwin13.4.0 (64-bit).
2. Semiparametric frailty models
The function survregbayes supports three commonly-used semiparametric frailty models:AFT, PH, and PO. The AFT model has survival and density functions S x ij ( t ) = S ( e x (cid:62) ij β + v i t ) , f x ij ( t ) = e x (cid:62) ij β + v i f ( e x (cid:62) ij β + v i t ) , (1) spBayesSurv, version 1.1.3 while the PH model has survival and density functions S x ij ( t ) = S ( t ) e x (cid:62) ij β + vi , f x ij ( t ) = e x (cid:62) ij β + v i S ( t ) e x (cid:62) ij β + vi − f ( t ) , (2)and the PO model has survival and density functions S x ij ( t ) = e − x (cid:62) ij β − v i S ( t )1 + ( e − x (cid:62) ij β − v i − S ( t ) , f x ij ( t ) = e − x (cid:62) ij β − v i f ( t )[1 + ( e − x (cid:62) ij β − v i − S ( t )] , (3)where β = ( β , . . . , β p ) (cid:62) is a vector of regression coefficients, v i is an unobserved frailtyassociated with s i , and S ( t ) is the baseline survival with density f ( t ) corresponding to x ij = and v i = 0. Let Γ( a, b ) denote a gamma distribution with mean a/b and N p ( µ , Σ )a p -variate normal distribution with mean µ and covariance Σ . The survregbayes functionimplements the following prior distributions: β ∼ N p ( β , S ) ,S ( · ) | α, θ ∼ TBP L ( α, S θ ( · )) , α ∼ Γ( a , b ) , θ ∼ N ( θ , V ) , ( v , . . . , v m ) (cid:62) | τ ∼ ICAR( τ ) , τ − ∼ Γ( a τ , b τ ) , or( v , . . . , v m ) (cid:62) | τ, φ ∼ GRF( τ , φ ) , τ − ∼ Γ( a τ , b τ ) , φ ∼ Γ( a φ , b φ ) , or( v , . . . , v m ) (cid:62) | τ ∼ IID( τ ) , τ − ∼ Γ( a τ , b τ )where TBP L , ICAR, GRF and IID refer to the transformed Bernstein polynomial (TBP)(Chen, Hanson, and Zhang 2014; Zhou and Hanson 2017) prior, intrinsic conditionally au-toregressive (ICAR) (Besag 1974) prior, Gaussian random field (GRF) prior, and independentGaussian (IID) prior distributions, respectively. The function argument prior allows usersto specify these prior parameters in a list with elements defined as follows:element maxL beta0 S0 a0 b0 theta0 V0 taua0 taub0 phia0 phib0 symbol L β S a b θ V a τ b τ a φ b φ We next briefly introduce these priors but leave details to Zhou and Hanson (2017).
TBP prior
In semiparametric survival analysis, a wide variety of Bayesian nonparametric priors can beused to model S ( · ); see M¨uller, Quintana, Jara, and Hanson (2015) and Zhou and Hanson(2015) for reviews. The TBP prior is attractive in that it is centered at a given paramet-ric family and it selects only smooth densities. For a fixed positive integer L , the priorTBP L ( α, S θ ( · )) is defined as S ( t ) = L (cid:88) j =1 w j I ( S θ ( t ) | j, L − j + 1) , w L ∼ Dirichlet( α, . . . , α ) , where w L = ( w , . . . , w L ) (cid:62) is a vector of positive weights, I ( ·| a, b ) denotes a beta cumulativedistribution function (cdf) with parameters ( a, b ), and { S θ ( · ) : θ ∈ Θ } is a parametricfamily of survival functions with support on positive reals R + . The log-logistic S θ ( t ) = { e θ t ) exp( θ ) } − , the log-normal S θ ( t ) = 1 − Φ { (log t + θ ) exp( θ ) } , and the Weibull S θ ( t ) = ournal of Statistical Software − exp (cid:8) − ( e θ t ) exp( θ ) (cid:9) families are implemented in survregbayes , where θ = ( θ , θ ) (cid:62) . Inour experience, the three centering distributions yield almost identical posterior inferencesbut in small samples one might be preferred. The random distribution S ( · ) is centered at S θ ( · ), i.e., E [ S ( t ) | α, θ ] = S θ ( t ). The parameter α controls how close the weights w j areto 1 /L , i.e., how close the shape of the baseline survival S ( · ) is relative to the prior guess S θ ( · ). Large values of α indicate a strong belief that S ( · ) is close to S θ ( · ); as α → ∞ , S ( · ) → S θ ( · ) with probability 1. Smaller values of α allow more pronounced deviations of S ( · ) from S θ ( · ). This adaptability makes the TBP prior attractive in its flexibility, but alsoanchors the random S ( · ) firmly about S θ ( · ): w j = 1 /L for j = 1 , . . . , L implies S ( t ) = S θ ( t )for t ≥
0. Moreover, unlike the mixture of Polya trees (Lavine 1992) or mixture of Dirichletprocess (Antoniak 1974) priors, the TBP prior selects smooth densities, leading to efficientposterior sampling.
ICAR and IID priors
For areal data, the ICAR prior smooths neighboring geographic-unit frailties v = ( v , . . . , v m ) (cid:62) .Let e ij be 1 if regions i and j share a common boundary and 0 otherwise; set e ii = 0. Thenthe m × m matrix E = [ e ij ] is called the adjacency matrix for the m regions. The priorICAR( τ ) on v is defined through the set of the conditional distributions v i |{ v j } j (cid:54) = i ∼ N m (cid:88) j =1 e ij v j /e i + , τ /e i + , i = 1 , . . . , m, (4)where e i + = (cid:80) mj =1 e ij is the number of neighbors of area s i . The induced prior on v underICAR is improper; the constraint (cid:80) mj =1 v j = 0 is used for identifiability (Banerjee, Carlin,and Gelfand 2014). Note that we assume that every region has at least one neighbor, so theproportionality constant for the improper density of v is ( τ − ) ( m − / (Lavine and Hodges2012).For non-spatial data, we consider the independent Gaussian prior IID( τ ), defined as v , v , . . . , v m iid ∼ N (0 , τ ) . (5) GRF priors
For georeferenced data, it is commonly assumed that v i = v ( s i ) arises from a Gaussian ran-dom field (GRF) { v ( s ) , s ∈ S} such that v = ( v , . . . , v m ) follows a multivariate Gaussiandistribution as v ∼ N m ( , τ R ), where τ measures the amount of spatial variation across lo-cations and the ( i, j ) element of R is modeled as R [ i, j ] = ρ ( s i , s j ). Here ρ ( · , · ) is a correlationfunction controlling the spatial dependence of v ( s ). In survregbayes the powered exponen-tial correlation function ρ ( s , s (cid:48) ) = ρ ( s , s (cid:48) ; φ ) = exp {− ( φ (cid:107) s − s (cid:48) (cid:107) ) ν } is used, where φ > ν ∈ (0 ,
2] is a pre-specifiedshape parameter which can be specified via prior$nu , and (cid:107) s − s (cid:48) (cid:107) refers to the distance(e.g., Euclidean, great-circle) between s and s (cid:48) . Therefore, the prior GRF( τ , φ ) is defined as v i |{ v j } j (cid:54) = i ∼ N − (cid:88) { j : j (cid:54) = i } p ij v j /p ii , τ /p ii , i = 1 , . . . , m, spBayesSurv, version 1.1.3 where p ij is the ( i, j ) element of R − . Full-scale approximation As m increases evaluating R − from R becomes computationally impractical. To overcomethis computational issue, we consider the FSA (Sang and Huang 2012) due to its capa-bility of capturing both large- and small-scale spatial dependence. Consider a fixed set of“knots” S ∗ = { s ∗ , . . . , s ∗ K } chosen from the study region. These knots are chosen usingthe function cover.design within the R package fields (Nychka, Furrer, Paige, and Sain2015), which computes space-filling coverage designs using the swapping algorithm (John-son, Moore, and Ylvisaker 1990). Let ρ ( s , s (cid:48) ) be the correlation between locations s and s (cid:48) .The usual predictive process approach (e.g., Banerjee et al. ρ ( s , s (cid:48) ) with ρ l ( s , s (cid:48) ) = ρ (cid:62) ( s , S ∗ ) ρ − KK ( S ∗ , S ∗ ) ρ ( s (cid:48) , S ∗ ), where ρ ( s , S ∗ ) = [ ρ ( s , s ∗ i )] Ki =1 is a K × ρ KK ( S ∗ , S ∗ ) = [ ρ ( s ∗ i , s ∗ j )] Ki,j =1 is a K × K correlation matrix at knots S ∗ . However, noting that ρ ( s , s (cid:48) ) = ρ l ( s , s (cid:48) ) + [ ρ ( s , s (cid:48) ) − ρ l ( s , s (cid:48) )], the predictive process discards entirely the residualpart ρ ( s , s (cid:48) ) − ρ l ( s , s (cid:48) ). In contrast, the FSA approach approximates the correlation function ρ ( s , s (cid:48) ) with ρ † ( s , s (cid:48) ) = ρ l ( s , s (cid:48) ) + ρ s ( s , s (cid:48) ) , (6)where ρ s ( s , s (cid:48) ) = { ρ ( s , s (cid:48) ) − ρ l ( s , s (cid:48) ) } ∆( s , s (cid:48) ) serves as a sparse approximate of the residualpart. Here ∆( s , s (cid:48) ) is a modulating function, which is specified so that ρ s ( s , s (cid:48) ) can well capturethe local residual spatial dependence while still permitting efficient computation. Motivatedby Konomi, Sang, and Mallick (2014), we first partition the total input space into B disjointblocks, and then specify ∆( s , s (cid:48) ) in a way such that the residuals are independent acrossinput blocks, but the original residual dependence structure within each block is retained.Specifically, the function ∆( s , s (cid:48) ) is taken to be 1 if s and s (cid:48) belong to the same block and 0otherwise. The approximated correlation function ρ † ( s , s (cid:48) ) in Equation 6 provides an exactrecovery of the true correlation within each block, and the approximation errors are ρ ( s , s (cid:48) ) − ρ l ( s , s (cid:48) ) for locations s and s (cid:48) in different blocks. Those errors are expected to be small formost entries because most of these location pairs are farther apart. To determine the blocks,we first use the R function cover.design to choose B ≤ m locations among the m locationsforming B blocks, then assign each s i to the block that is closest to s i . Here B does notneed to be equal to K . When B = 1, no approximation is applied to the correlation ρ .When B = m , it reduces to the approach of Finley et al. (2009), so the local residual spatialdependence may not be well captured.Applying the above FSA approach to approximate the correlation function ρ ( s , s (cid:48) ), we canapproximate the correlation matrix R with ρ † mm = ρ l + ρ s = ρ mK ρ − KK ρ (cid:62) mK + (cid:16) ρ mm − ρ mK ρ − KK ρ (cid:62) mK (cid:17) ◦ ∆ , (7)where ρ mK = [ ρ ( s i , s ∗ j )] i =1: m,j =1: K , ρ KK = [ ρ ( s ∗ i , s ∗ j )] Ki,j =1 , and ∆ = [∆( s i , s j )] mi,j =1 . Here, thenotation “ ◦ ” represents the element-wise matrix multiplication. To avoid numerical instability,we add a small nugget effect (cid:15) = 10 − when defining R , that is, R = (1 − (cid:15) ) ρ mm + (cid:15) I m . Itfollows from Equation 7 that R can be approximated by R † = (1 − (cid:15) ) ρ † mm + (cid:15) I m = (1 − (cid:15) ) ρ mK ρ − KK ρ (cid:62) mK + R s , ournal of Statistical Software R s = (1 − (cid:15) ) (cid:0) ρ mm − ρ mK ρ − KK ρ (cid:62) mK (cid:1) ◦ ∆ + (cid:15) I m . Applying the Sherman-Woodbury-Morrison formula for inverse matrices, we can approximate R − by (cid:16) R † (cid:17) − = R − s − (1 − (cid:15) ) R − s ρ mK (cid:104) ρ KK + (1 − (cid:15) ) ρ (cid:62) mK R − s ρ mK (cid:105) − ρ (cid:62) mK R − s . (8)In addition, the determinant of R can be approximated bydet (cid:16) R † (cid:17) = det (cid:110) ρ KK + (1 − (cid:15) ) ρ (cid:62) mK R − s ρ mK (cid:111) det( ρ KK ) − det( R s ) . (9)Since the m × m matrix R s is a block matrix, the right-hand sides of Equations 8 and 9involve only inverses and determinants of K × K low-rank matrices and m × m block diagonalmatrices. Thus the computational complexity can be greatly reduced relative to the expensivecomputational cost of using original correlation function for large value of m . However, forsmall m , e.g., m < R due to thecomplexity of FSA’s implementation. Note that K and B can be specified via prior$K and prior$B , respectively. The likelihood function for ( w L , θ , β , v ) is given by L ( w L , θ , β , v ) = m (cid:89) i =1 n i (cid:89) j =1 (cid:2) S x ij ( a ij ) − S x ij ( b ij ) (cid:3) I { a ij
043 patients (Henderson et al.
LeukSurv in the package. It is of interest to investigate possiblespatial variation in survival after accounting for known subject-specific prognostic factors,which include age , sex , white blood cell count ( wbc ) at diagnosis, and the Townsend score( tpi ) for which higher values indicates less affluent areas. Both exact residential locations ofall patients and their administrative districts (the boundary file is named as nwengland.bnd in the package) are available, so we can fit both geostatistical and areal models. PO model with ICAR frailties
If the IID or ICAR frailties are considered, to easily identify the correspondence betweenfrailties and clusters/regions, we program the function survregbayes so that the input datasetshould be sorted by the cluster variable before any use. The following code is used to sortthe dataset by district and obtain the adjacency matrix E . R> library("coda")R> library("survival")R> library("spBayesSurv")R> library("fields")R> library("BayesX")R> library("R2BayesX")R> data("LeukSurv")R> d <- LeukSurv[order(LeukSurv$district), ]R> head(d) time cens xcoord ycoord age sex wbc tpi district24 1 1 0.4123484 0.4233738 44 1 281.0 4.87 162 3 1 0.3925028 0.4531422 72 1 0.0 7.10 1 ournal of Statistical Software
68 4 1 0.4167585 0.4520397 68 0 0.0 5.12 1128 9 1 0.4244763 0.4123484 61 1 0.0 2.90 1129 9 1 0.4145535 0.4520397 26 1 0.0 6.72 1163 15 1 0.4013230 0.4785006 67 1 27.9 1.50 1
R> nwengland <- read.bnd(system.file("otherdata/nwengland.bnd",+ package = "spBayesSurv"))R> adj.mat <- bnd2gra(nwengland)R> E <- diag(diag(adj.mat)) - as.matrix(adj.mat)
The following code is used to fit the PO model with ICAR frailties using the TBP prior with L = 15 and default settings for other priors. A burn-in period of 5,000 iterates was consideredand the Markov chain was subsampled every 5 iterates to get a final chain size of 2,000. Theargument ndisplay = 1000 will display the number of saved scans after every 1,000 savediterates. If the argument InitParamMCMC = TRUE (not used here as it is the default setting),then an initial chain with nburn = 5000 , nsave = 5000 , nkip = 0 and ndisplay = 1000 will be run under parametric models; otherwise, the initial values are obtained from fittingparametric non-frailty models via survreg . The total running time is 166 seconds. R> set.seed(1)R> mcmc <- list(nburn = 5000, nsave = 2000, nskip = 4, ndisplay = 1000)R> prior <- list(maxL = 15)R> ptm <- proc.time()R> res1 <- survregbayes(formula = Surv(time, cens) ~ age + sex + wbc + tpi ++ frailtyprior("car", district), data = d, survmodel = "PO",+ dist = "loglogistic", mcmc = mcmc, prior = prior, Proximity = E)R> proc.time() - ptm user system elapsed165.919 0.296 166.354
The term frailtyprior("car", district) indicates that the ICAR prior in Equation 4is used. One can also incorporate the IID prior in Equation 5 via frailtyprior("iid",district) . The non-frailty model can be fit by removing the frailtyprior term. Theargument survmodel is used to indicate which model will be fit; choices include "PH" , "PO" ,and "AFT" . The argument dist is used to specify the distribution family of S θ ( · ) definedin Section 2.1, and the choices include "loglogistic" , "lognormal" , and "weibull" . Theargument prior is used to specify user-defined hyperparameters, e.g., for p = 3, L = 15, β = , S = 10 I p , θ = , V = 10 I , a = b = 1, and a τ = b τ = 1, the prior can bespecified as below. R> prior <- list(maxL = 15, beta0 = rep(0, 3), S0 = diag(10, 3),+ theta0 = rep(0, 2), V0 = diag(10, 2), a0 = 1, b0 = 1,+ taua0 = 1, taub0 = 1) If prior = NULL , then the default hyperparameters given in Section 2.2 would be used. Noteby default survregbayes standardizes each covariate by subtracting the sample mean and0 spBayesSurv, version 1.1.3 dividing the sample standard deviation. Therefore, the user-specified hyperparameters shouldbe based on the model with scaled covariates unless the argument scale.designX = FALSE is added.The output from applying the summary function to the returned object res1 is given below. R> (sfit1 <- summary(res1))
Proportional Odds model:Call:survregbayes(formula = Surv(time, cens) ~ age + sex + wbc + tpi +frailtyprior("car", district), data = d, survmodel = "PO",dist = "loglogistic", mcmc = mcmc, prior = prior, Proximity = E)Posterior inference of regression coefficients(Adaptive M-H acceptance rate: 0.2731):Mean Median Std. Dev. 95%CI-Low 95%CI-Uppage 0.0519835 0.0518955 0.0034329 0.0455544 0.0589767sex 0.1238558 0.1241657 0.1061961 -0.0854203 0.3274537wbc 0.0059439 0.0059223 0.0008163 0.0043996 0.0074789tpi 0.0598826 0.0597254 0.0159244 0.0286519 0.0904957Posterior inference of conditional CAR frailty varianceMean Median Std. Dev. 95%CI-Low 95%CI-Uppvariance 0.080346 0.056350 0.082950 0.001709 0.299395Log pseudo marginal likelihood: LPML=-5925.194Deviance Information Criterion: DIC=11849.82Watanabe-Akaike information criterion: WAIC=11850.39Number of subjects: n=1043
We can see that age , wbc and tpi are significant risk factors for leukemia survival. Forexample, lower age decreases the odds of a patient dying by any time; holding other predictorsconstant, a 10-year decrease in age cuts the odds of dying by exp( − × . ≈ τ is 0 .
08. The LPML, DIC and WAIC are -5925, 11850 and 11850,respectively.The following code is used to produce trace plots (Figure 1) for β and τ . Note that themixing for τ is not very satisfactory. This is not surprising, since we are using very vaguegamma prior Γ(0 . , . ,
1) on τ or run a longer chain with higher thin to improvethe mixing. R> par(mfrow = c(3, 2))R> par(cex = 1, mar = c(2.5, 4.1, 1, 1))R> traceplot(mcmc(res1$beta[1,]), xlab = "", main = "age")R> traceplot(mcmc(res1$beta[2,]), xlab = "", main = "sex")R> traceplot(mcmc(res1$beta[3,]), xlab = "", main = "wbc") ournal of Statistical Software . . age − . . . . sex . . . wbc . . . tpi . . . . . tau^2 Figure 1: Leukemia survival data. Trace plots for β , τ and α under the PO model withICAR frailties. R> traceplot(mcmc(res1$beta[4,]), xlab = "", main = "tpi")R> traceplot(mcmc(res1$tau2), xlab = "", main = "tau^2")
The code below is used to generate the Cox-Snell plots with 10 posterior residuals (Figure 2,panel a).
R> set.seed(1)R> cox.snell.survregbayes(res1, ncurves = 10)
The code below is used to generate survival curves for female patients with wbc =38.59 and tpi =0.3398 at different ages (Figure 2, panel b).
R> tgrid <- seq(0.1, 5000, length.out = 300);R> xpred <- data.frame(age = c(49, 65, 74), sex = c(0, 0, 0),+ wbc = c(38.59, 38.59, 38.59), tpi = c(0.3398, 0.3398, 0.3398),+ row.names = c("age=49", "age=65", "age=74"))R> plot(res1, xnewdata = xpred, tgrid = tgrid, cex = 2)
The code below is used to generate the map of posterior means of frailties for each district(Figure 2, panel c). Note that the posterior median of frailties can be extracted similarly byreplacing mean below with median in the apply function.
R> frail0 <- apply(res1$v, 1, mean)R> frail <- frail0[as.integer(names(nwengland))]R> values <- cbind(as.integer(names(nwengland)), frail)R> op <- par(no.readonly = TRUE)R> par(mar = c(3, 0, 0, 0))R> plotmap(nwengland, x = values, col = (gray.colors(10, 0.3, 1))[10:1],+ pos = "bottomleft", width = 0.5, height = 0.04) spBayesSurv, version 1.1.3 (a) . . . . . . time s u r v i v a l age=49age=65age=74 (b) −0.23 0 0.23 (c) Figure 2: Leukemia survival data. PO model with ICAR frailties. (a) Cox-Snell plot. (b)Survival curves with 95% credible interval bands for female patients with wbc =38.59 and tpi =0.3398 at different ages. (c) Map for the posterior mean frailties; larger frailties meanhigher mortality rate overall.
PO model with GRF frailties
Note that all coordinates are distinct, so we have m = 1043 and n i = 1 in terms of ournotation. To use frailtyprior to specify the prior, we need to create an ID variable consistingof 1043 distinct values. The powered exponential correlation function with ν = 1 is used. Tospecify the number of knots and blocks for the FSA of R , we consider K = 100 and B = 1043.The code below is used to fit a PO model with GRF frailties under above settings. The runningtime is a bit under three hours. R> set.seed(1)R> mcmc <- list(nburn = 5000, nsave = 2000, nskip = 4, ndisplay = 1000)R> prior <- list(maxL = 15, nu = 1, nknots = 100, nblock = 1043)R> d$ID <- 1:nrow(d)R> locations <- cbind(d$xcoord, d$ycoord);R> ptm <- proc.time()R> res2 <- survregbayes(formula = Surv(time, cens) ~ age + sex + wbc + tpi ++ frailtyprior("grf", ID), data = d, survmodel = "PO",+ dist = "loglogistic", mcmc = mcmc, prior = prior,+ Coordinates = locations)R> proc.time() - ptm user system elapsed10079.006 97.039 10176.650
R> (sfit2 <- summary(res2))
Posterior inference of regression coefficients(Adaptive M-H acceptance rate: 0.2726):Mean Median Std. Dev. 95%CI-Low 95%CI-Uppage 0.0526668 0.0527180 0.0034351 0.0460261 0.0596917 ournal of Statistical Software sex 0.1310119 0.1318825 0.1069728 -0.0748948 0.3457847wbc 0.0060590 0.0060293 0.0008156 0.0044876 0.0077388tpi 0.0606026 0.0609221 0.0158076 0.0300292 0.0918792Posterior inference of frailty varianceMean Median Std. Dev. 95%CI-Low 95%CI-Uppvariance 0.06179 0.05290 0.03261 0.02376 0.14086Posterior inference of correlation function range phiMean Median Std. Dev. 95%CI-Low 95%CI-Upprange 19.138 17.245 7.305 8.701 35.094Log pseudo marginal likelihood: LPML=-5923.402Deviance Information Criterion: DIC=11845.78Watanabe-Akaike information criterion: WAIC=11846.78Number of subjects: n=1043 The trace plots for β , τ and φ (Figure 3), Cox-Snell residuals and survival curves (Figure 4)can be obtained using the same code used for the PO model with ICAR frailties. The codebelow is used to generate the map of posterior means of frailties for each location (Figure 4). R> frail <- round(apply(res2$v, 1, mean), 3)R> nclust <- 5R> frail.cluster <- cut(frail, breaks = nclust)R> frail.names <- names(table(frail.cluster))R> rbPal <- colorRampPalette(c( ' blue ' , ' red ' ))R> frail.colors <- rbPal(nclust)[as.numeric(frail.cluster)]R> par(mar = c(3, 0, 0, 0))R> plot(nwengland)R> points(cbind(d$xcoord,d$ycoord), col = frail.colors)R> legend("topright", title = "frailty values", legend = frail.names,+ col = rbPal(nclust), pch = 20, cex = 1.7) Note that the mixing for τ and φ is very poor. This may be partly due to the fact we areupdating large dimensional ( m = 1 , Let x = ( x , . . . , x p ) (cid:62) denote the p -vector of covariates in general. The most direct approachis to multiply β (cid:96) by a latent Bernoulli variable γ (cid:96) for (cid:96) = 1 , . . . , p , where γ (cid:96) = 1 indicates thepresence of covariate x (cid:96) in the model, and then assume an appropriate prior on ( β , γ ), where γ = ( γ , . . . , γ p ) (cid:62) . Following Kuo and Mallick (1998) and Hanson, Branscum, and Johnson4 spBayesSurv, version 1.1.3 . . . age − . . . . sex . . . wbc . . . tpi . . tau^2 phi Figure 3: Leukemia survival data. Trace plots for β , τ and α under the PO model with GRFfrailties. (a) . . . . . . time s u r v i v a l age=49age=65age=74 (b) lllllllllllllllll llll llllllll lllllllllllllllllllllll llllllllllll lllllllllllllllllllllll llllllllllllllll l lllll lllllllllllllllllll lll lll llll llll lllllll ll llll lll llllll ll lllllll llllll llll llll llllllllll llllllll ll ll lllllllll lllll lll ll l lllll lllll llll llllll llll llllll ll llllllll lll ll llll l ll llll llll ll l lllllll llll lllll lllllll lll l llll lllll ll ll llllllllllll llllll l ll llll ll llll lllllll ll lllllll ll ll llll llll lll llllll lll ll lllll ll lll lllllll ll ll ll lllllll llll llllll lllll lllllllllllll lll llll lll l ll ll ll l ll l lll ll ll l ll ll lll ll l l lllllll lll l llllll llllllllllll lllll llllll llllllllllllllllllllllll lllll lll lllll lll lll llllll lllll lll llll lllll ll llll l ll lll ll ll lll lllll ll llllll llll lll llllllll llll llll ll llllllllll ll lllllll llll llll lllll ll llll ll llllll llllllll llll llllll lll lll ll llllllllllll lll ll ll llllll lll ll l llll l ll lll ll llll llll llll lll llll ll lll lllll ll lll ll lllll ll lllllllllllllllllllll lllllll lllllll ll lll lll lll ll ll lllllll ll lllllllll l lll lll llll l llllll llll llll lllllllllllll lll llllll ll llll lll l llllll l llllll ll llll lll l lllllll ll l ll lllllll l lllll llll lll ll ll llll l llll ll ll lll ll lll l ll lllllllll l lllll frailty values(−0.32,−0.208](−0.208,−0.0978](−0.0978,0.0128](0.0128,0.123](0.123,0.235] (c) Figure 4: Leukemia survival data. PO model with GRF frailties. (a) Cox-Snell plot. (b)Survival curves with 95% credible interval bands for female patients with wbc =38.59 and tpi =0.3398 at different ages. (c) Map for the posterior mean frailties; larger frailties meanhigher mortality rate overall. ournal of Statistical Software γ , . . . , γ p iid ∼ Bern(0 .
5) and β ∼ N p ( , gn ( X (cid:62) X ) − ) , where X is the usual design matrix, but with mean-centered covariates, i.e., (cid:62) n X = (cid:62) p , and g is chosen by picking a number M such that a random e x (cid:62) β is less than M with probability q , i.e., approximately g = (cid:2) log M / Φ − ( q ) (cid:3) /p . The function survregbayes sets M = 10 and q = 0 . M and q via prior$M and prior$q ,respectively. The MCMC procedure is described in Zhou and Hanson (2017).To perform variable selection for the leukemia survival data, we simply need to add theargument selection=TRUE to the function survregbayes . A part of the output from summary is also shown. The model with age , wbc and tpi has the highest proportion (89.8%), andthus can be served as the final model. R> set.seed(1)R> res3 <- survregbayes(formula = Surv(time, cens) ~ age + sex + wbc + tpi ++ frailtyprior("car", district), data = d, survmodel = "PO",+ dist = "loglogistic", mcmc = mcmc, prior = prior, Proximity = E,+ selection = TRUE)R> (sfit3 <- summary(res3))
Variable selection:age,wbc,tpi age,sex,wbc,tpi age,wbcprop. 0.8975 0.1010 0.0015
Many authors have found parametric models to fit as well or better than competing semipara-metric models (Cox and Oakes 1984, p. 123; Nardi and Schemper 2003). The semiparametric– or more accurately richly parametric – formulation of the AFT, PH and PO models pre-sented here have their baseline survival functions centered at a parametric family S θ ( t ). Notethat z J − = implies S ( t ) = S θ ( t ). Therefore, testing H : z J − = versus H : z J − (cid:54) = leads to the comparison of the semiparametric model with the underlying parametric model.Let BF be the Bayes factor between H and H . Zhou et al. (2017) proposed to esti-mate BF by a large-sample approximation to the generalized Savage-Dickey density ratio(Verdinelli and Wasserman 1995). Adapting their approach BF is estimated (cid:100) BF = p ( | ˆ α ) N J − ( ; ˆ m , ˆ Σ ) , where p ( | α ) = Γ( αJ ) / [ J α Γ( α )] J is the prior density of z J − evaluated at z J − = , ˆ α isthe posterior mean of α , N p ( · ; m , Σ ) denotes a p -variable normal density with mean m andcovariance Σ , and ˆ m and ˆ Σ are posterior mean and covariance of z J − .The Bayes factor BF under the semiparametric PO model with ICAR frailties can beobtained using the code below (here the object res1 is obtained in Section 2.4). R> BF.survregbayes(res1)[1] 82.12799 spBayesSurv, version 1.1.3 The BF = 82 > survregbayes also supports the efficient fitting of parametric frailty modelswith loglogistic, lognormal or Weibull baseline functions. In parametric models, the prior for θ can be set to be relatively vague. Setting a at any negative value will force the α to befixed at the value specified in the argument state . For example, setting prior <- list(a0= -1) and state = list(alpha = 1) will fix α = 1 throughout the MCMC; setting prior =list(a0 = -1) and state = list(alpha = Inf) will fit a parametric model. The followingcode fits a parametric loglogistic PO model with ICAR frailties to the leukemia survival data.The LPML is -5950, much worse than the value under the semiparametric PO model. R> set.seed(1)R> prior <- list(maxL = 15, a0 = -1, thete0 = rep(0, 2), V0 = diag(1e10, 2))R> state <- list(alpha = Inf)R> ptm <- proc.time()R> res11 <- survregbayes(formula = Surv(time, cens) ~ age + sex + wbc + tpi ++ frailtyprior("car", district), data = d, survmodel = "PO",+ dist = "loglogistic", mcmc = mcmc, prior = prior, state = state,+ Proximity = E, InitParamMCMC = FALSE)R> proc.time() - ptm user system elapsed25.037 0.115 25.239
R> (sfit11 <- summary(res11))
Proportional Odds model:Call:survregbayes(formula = Surv(time, cens) ~ age + sex + wbc + tpi +frailtyprior("car", district), data = d, survmodel = "PO",dist = "loglogistic", mcmc = mcmc, prior = prior, state = state,Proximity = E, InitParamMCMC = FALSE)Posterior inference of regression coefficients(Adaptive M-H acceptance rate: 0.2844):Mean Median Std. Dev. 95%CI-Low 95%CI-Uppage 0.0504253 0.0504362 0.0033318 0.0439945 0.0568477sex 0.1187297 0.1134544 0.1109109 -0.0912841 0.3374972wbc 0.0062192 0.0062147 0.0007395 0.0048068 0.0076600tpi 0.0602207 0.0603376 0.0156038 0.0299010 0.0915584Posterior inference of conditional CAR frailty varianceMean Median Std. Dev. 95%CI-Low 95%CI-Uppvariance 0.078627 0.055078 0.082164 0.002005 0.305202Log pseudo marginal likelihood: LPML=-5949.919Deviance Information Criterion: DIC=11899.46 ournal of Statistical Software Watanabe-Akaike information criterion: WAIC=11899.84Number of subjects: n=1043
The survival time t ij is left-truncated at u ij ≥ u ij is the time when the ij th subject isfirst observed. Left-truncation often occurs when age is used as the time scale. Given theobserved left-truncated data { ( u ij , a ij , b ij , x ij , s i ) } , where a ij ≥ u ij , the likelihood function inEquation 10 becomes L ( w J , θ , β , v ) = m (cid:89) i =1 n i (cid:89) j =1 (cid:2) S x ij ( a ij ) − S x ij ( b ij ) (cid:3) I { a ij a ij ) = P ( t ij > a ij | t ij > t ij,o ij ) o ij − (cid:89) k =1 P ( t ij > t ij,k +1 | t ij > t ij,k )= S x ij,oij ( a ij ) S x ij,oij ( t ij,o ij ) o ij − (cid:89) k =1 S x ij,k ( t ij,k +1 ) S x ij,k ( t ij,k ) . Thus one can replace the observation ( u ij , a ij , b ij , x ij ( t ) , s i ) by a set of new o ij observations( t ij, , t ij, , ∞ , x ij, , s i ), ( t ij, , t ij, , ∞ , x ij, , s i ), . . . , ( t ij,o ij , a ij , b ij , x ij,o ij , s i ). This way we get anew left-truncated data set of size (cid:80) mi =1 (cid:80) n i j =1 o ij . Then the likelihood function becomes L ( w J , θ , β , v ) = m (cid:89) i =1 n i (cid:89) j =1 (cid:26) (cid:104) S x ij,oij ( a ij ) − S x ij,oij ( b ij ) (cid:105) I { a ij
We use the primary biliary cirrhosis (PBC) dataset (available in the package survival as pbc ) as an example to show how to incorporate time-dependent covariates in the function8 spBayesSurv, version 1.1.3 survregbayes . Although this is not a spatial dataset, spatial frailties can be added similarlyas in Section 2.4. The following code is copied from Therneau, Crowson, and Atkinson (2017)to create the data frame with time-dependent covariates. R> temp <- subset(pbc, id <= 312, select = c(id:sex, stage)) id tstart tstop endpt bili protime1 1 0 192 0 14.5 12.22 1 192 400 2 21.3 11.23 2 0 182 0 1.1 10.64 2 182 365 0 0.8 11.05 2 365 768 0 1.0 11.66 2 768 1790 0 1.9 10.6
We can fit the Bayesian PH model with TBP baseline as follows. The output for regressioncoefficients is partial.
R> set.seed(1)R> mcmc <- list(nburn = 5000, nsave = 2000, nskip = 4, ndisplay = 1000)R> ptm <- proc.time()R> fit1 <- survregbayes(Surv(tstart, tstop, endpt == 2) ~ log(bili) ++ log(protime), data = pbc2, survmodel = "PH", dist = "loglogistic",+ mcmc = mcmc, subject.num = id)R> proc.time() - ptm user system elapsed227.626 0.434 228.243
R> summary(fit1)
Proportional hazards model:Call:survregbayes(formula = Surv(tstart, tstop, endpt == 2) ~ log(bili) +log(protime), data = pbc2, survmodel = "PH", dist = "loglogistic",mcmc = mcmc, subject.num = id)Posterior inference of regression coefficients(Adaptive M-H acceptance rate: 0.2135):Mean Median Std. Dev. 95%CI-Low 95%CI-Upplog(bili) 1.29937 1.30058 0.09452 1.11354 1.48584log(protime) 4.18500 4.20421 0.37052 3.43850 4.84161 ournal of Statistical Software Log pseudo marginal likelihood: LPML=-1018.001Deviance Information Criterion: DIC=2032.754Watanabe-Akaike information criterion: WAIC=2035.861Number of subjects: n=1807
Equivalently, one can also run the following code to obtain the same analysis. The argument truncation_time is used to specify the start time point for each time interval, i.e., tstart .The end time point tstop together with endpt are formulated as interval censored data using type = "interval2" of Surv . This format is more general than the former one, as one caneasily incorporate interval censored data.
R> pbc2$tleft <- pbc2$tstop; pbc2$tright <- pbc2$tstop;R> pbc2$tright[which(pbc2$endpt! = 2)] <- NA;R> fit11 <- survregbayes(Surv(tleft, tright, type = "interval2") ~ log(bili) ++ log(protime), data = pbc2, survmodel = "PH", dist = "loglogistic",+ mcmc = mcmc, truncation_time = tstart, subject.num = id);
3. GAFT frailty models
The generalized accelerated failure time (GAFT) frailty model (Zhou et al. S ( t ) to depend oncertain covariates, say a q -dimensional vector z ij which is usually a subset of x ij . Specifically,the GAFT frailty model is given by S x ij ( t ) = S , z ij (cid:16) e − x (cid:62) ij β − v i t (cid:17) , or equivalently, y ij = log( t ij ) = ˜ x (cid:62) ij ˜ β + v i + (cid:15) ij , where ˜ x ij = (1 , x (cid:62) ij ) (cid:62) includes an intercept, ˜ β = ( β , β (cid:62) ) (cid:62) is a vector of corresponding coeffi-cients, (cid:15) ij is a heteroscedastic error term independent of v i , and P ( e β + (cid:15) ij > t | z ij ) = S , z ij ( t ).Note the regression coefficients β here are defined differently with those in Equation 1. Herewe assume (cid:15) ij | G z ij ind. ∼ G z ij , where G z is a probability measure defined on R for every z ∈ X ; this defines a model for theentire collection of probability measures G X = { G z : z ∈ X } so that each element is allowedto smoothly change with the covariates z . The frailtyGAFT function considers the followingprior distributions: ˜ β ∼ N p +1 ( m , S ) G z | α, σ ∼ LDTFP L ( α, σ ) , α ∼ Γ( a , b ) , σ − ∼ Γ( a σ , b σ ) , ( v , . . . , v m ) (cid:62) | τ ∼ ICAR( τ ) , τ − ∼ Γ( a τ , b τ ) , or( v , . . . , v m ) (cid:62) | τ, φ ∼ GRF( τ , φ ) , τ − ∼ Γ( a τ , b τ ) , φ ∼ Γ( a φ , b φ ) , or( v , . . . , v m ) (cid:62) | τ ∼ IID( τ ) , τ − ∼ Γ( a τ , b τ )0 spBayesSurv, version 1.1.3 where LDTFP L refers to the linear dependent tailfree process prior (LDTFP) prior as de-scribed in (Zhou et al. prior allows users to specify theseprior parameters in a list with elements defined as follows:element maxL m0 S0 a0 b0 siga0 sigb0 taua0 taub0 phia0 phib0 symbol L m S a b a σ b σ a τ b τ a φ b φ The LDTFP prior considered in Zhou et al. (2017) is centered at a normal distribution Φ σ with mean 0 and variance σ , that is, E ( G z ) = Φ σ for every z ∈ X . Define the function k σ ( x ) = (cid:100) L Φ σ ( x ) (cid:101) , where (cid:100) x (cid:101) is the ceiling function, the smallest integer greater than orequal to x . Further define probability p z ( k ) for k = 1 , . . . , L as p z ( k ) = L (cid:89) l =1 Y l, (cid:100) k l − L (cid:101) ( z ) , where Y j +1 , k − ( z ) = (cid:0) {− ˜ z (cid:62) γ j,k } (cid:1) − and Y j +1 , k ( z ) = 1 − Y j +1 , k − ( z ) for j =0 , . . . , L − k = 1 , . . . , j , where ˜ z = (1 , z (cid:62) ) (cid:62) includes an intercept, and γ j,k = ( γ j,k, , . . . , γ j,k,q ) (cid:62) is a vector of coefficients. Note there are 2 L − γ = { γ j,k } , e.g.,for L = 3, γ = { γ , , γ , , γ , , γ , , γ , , γ , , γ , } . For a fixed integer L >
0, the randomdensity associated with LDTFP L ( α, σ ) is defined as f z ( e ) = 2 L φ σ ( e ) p z { k σ ( e ) } , γ j,k ind. ∼ N q +1 (cid:18) , nα ( j + 1) ( Z (cid:62) Z ) − (cid:19) with cdf G z ( e ) = p z { k σ ( e ) } (cid:8) L Φ σ ( e ) − k σ ( e ) (cid:9) + k σ ( e ) (cid:88) k =1 p z ( k ) , (11)where Z is the n × ( q + 1) design matrix with mean-centered covariates ˜ z ij s. Furthermore,the LDTFP is specified by setting γ , ≡ , such that for every z ∈ X , G z is almost surely amedian-zero probability measure.The function frailtyGAFT sets the following hyperparameters as defaults: m = , S =10 I p +1 , a = b = 1, a τ = b τ = 1, and a σ = 2 + ˆ σ / (100ˆ v ), b σ = ˆ σ ( a σ − σ and ˆ v are the estimates of σ and its asymptotic variance from fitting the parametriclognormal AFT model, respectively. Note here we assume a somewhat informative prior on σ so that its mean is ˆ σ and variance is 100ˆ v . For the GRF prior, we again set a φ = 2 and b φ = ( a φ − /φ so that the prior of φ has mode at φ and the prior mean of 1 /φ is 1 /φ withinfinite variance. Here φ satisfies ρ ( s (cid:48) , s (cid:48)(cid:48) ; φ ) = 0 . (cid:107) s (cid:48) − s (cid:48)(cid:48) (cid:107) = max ij (cid:107) s i − s j (cid:107) .Note by default frailtyGAFT standardizes each covariate by subtracting the sample meanand dividing the sample standard deviation. Therefore, the user-specified hyperparametersshould be based on the model with scaled covariates unless the argument scale.designX =FALSE is added. The GAFT frailty model includes the following as important special cases: an AFT frailtymodel with nonparametric baseline where G z = G z (cid:48) for all z = z (cid:48) and parametric baseline ournal of Statistical Software G z = Φ σ for all z ∈ X . Hypothesis tests can be constructed based on the LDTFPcoefficients { γ l,k : k = 1 , . . . , l , l = 1 , . . . , L − } , where γ l,k = ( γ l,k, , . . . , γ l,k,q ) (cid:62) . Let γ l,k, − j denote the subvector of γ l,k without element γ l,k,j for j = 0 , . . . , q . Set Υ j = ( γ l,k,j , k =1 , . . . , l , l = 1 , . . . , L − (cid:62) , Υ − j = ( γ (cid:62) l,k, − j , k = 1 , . . . , l , l = 1 , . . . , L − (cid:62) and Υ = ( γ (cid:62) l,k , k =1 , . . . , l , l = 1 , . . . , L − (cid:62) . Testing the hypotheses H : Υ − = and H : Υ = leadsto global comparisons of the proposed model with the above two special cases respectively.Similarly, we may also test the null hypothesis H : Υ j = for the j th covariate effect of z on the baseline survival, j = 1 , . . . , q .Suppose we wish to test H : Υ j = versus H : Υ j (cid:54) = , for fixed j ∈ { , . . . , q } . FollowingZhou et al. (2017), the Bayes factor between hypotheses H and H can be approximated byˆ BF = L − (cid:89) l =1 2 l (cid:89) k =1 N (cid:18) (cid:12)(cid:12)(cid:12)(cid:12) , n ˆ α ( l + 1) ( Z (cid:62) Z ) − jj (cid:19) N L − ( Υ j = ; ˆ m j , ˆ S j ) , where N p ( · ; m , S ) denotes a p -variate normal density with mean m and covariance matrix S ,and ˆ m j and ˆ S j are the sample mean and covariance for Υ j . The code below is used to fit the GAFT model with ICAR frailties for the leukemia survivaldata. As suggested by Zhou et al. (2017), the gamma prior Γ( a = 5 , b = 1) is used for α .We include all four covariates in modeling the baseline survival function. R> set.seed(1)R> mcmc <- list(nburn = 5000, nsave = 2000, nskip = 4, ndisplay = 1000)R> prior <- list(maxL = 4, a0 = 5, b0 = 1)R> ptm <- proc.time()R> res1 <- frailtyGAFT(formula = Surv(time, cens) ~ age + sex + wbc + tpi ++ baseline(age, sex, wbc, tpi) + frailtyprior("car", district),+ data = d, mcmc = mcmc, prior = prior, Proximity = E)R> (sfit1 <- summary(res1))
Generalized accelerated failure time frailty model:Call:frailtyGAFT(formula = Surv(time, cens) ~ age + sex + wbc + tpi +baseline(age, sex, wbc, tpi) + frailtyprior("car", district),data = d, mcmc = mcmc, prior = prior, Proximity = E)Posterior inference of regression coefficientsMean Median Std. Dev. 95%HPD-Low 95%HPD-Uppintercept 8.589761 8.607783 0.288535 7.982265 9.140489age -0.051342 -0.051508 0.003987 -0.058561 -0.041985sex -0.267978 -0.288883 0.164310 -0.533180 0.064909wbc -0.004161 -0.004322 0.001001 -0.005931 -0.001864 spBayesSurv, version 1.1.3 . . . intercept − . − . − . age − . − . − . . . sex − . − . − . wbc − . − . − . tpi . . . . . . tau^2 (a) . . . . . . time s u r v i v a l age=49age=65age=74 (b) −0.52 0 0.52 (c) Figure 5: Leukemia survival data. GAFT model with ICAR frailties. (a) Trace plots for β , τ and α . (b) Survival curves with 95% credible interval bands for female patients with wbc =38.59 and tpi =0.3398 at different ages. (c) Map for the negative posterior mean frailties;larger values mean higher mortality rate overall. tpi -0.065335 -0.067061 0.019601 -0.099992 -0.023739Bayes factors for LDTFP covariate effects:intercept age sex wbc tpi overall normality220.2500 16.2494 1.1579 28.1776 0.4842 11.2454 1787.3269Log pseudo marginal likelihood: LPML=-5937.304Number of subjects:=1043 R> proc.time() - ptm user system elapsed444.393 2.433 454.270
The Bayes factors for testing age and wbc effects on LDTFP are 16 and 28, respectively,indicating that the baseline survival function under the AFT model depends on age and wbc , and thus GAFT should be considered. The trace plots, survival curves and frailty map(Figure 5) can be obtained using the code similarly as in Section 2.4. The only differencefor plotting survival curves is that we need to specify the baseline covariates by including theargument xtfnewdata = xpred into the plot function. Note that the mixing for covariateeffects is okay but not great due to the non-smoothness of Polya trees. In this case, we needto run a longer chain with much higher thinning as suggested in Zhou et al. (2017).
4. Survival models via spatial copulas
In environmental studies, survival times (e.g. time to water pollution) often present a strongspatial dependence after adjusting for available risk factors, making frailty models extremelydifficult to fit because of the strong posterior dependency among frailties. The spatial copula ournal of Statistical Software n i = 1), right-censored spatial data. Suppose subjects are observed at n distinct spa-tial locations s , . . . , s n . Let t i be a random event time associated with the subject at s i and x i be a related p -dimensional vector of covariates, i = 1 , . . . , n . For right-censored data, we onlyobserve t oi and a censoring indicator δ i for each subject, where δ i equals 1 if t oi = t i and equals0 if t i is censored at t oi . Therefore, the observed data will be D = { ( t oi , δ i , x i , s i ); i = 1 , . . . , n } .Note although the models below are developed for spatial survival data, non-spatial data arealso accommodated.In the context of survival models, the idea of spatial copula approach is to first assume thatthe survival time t i at location s i marginally follows a model S x i ( t ), then model the jointdistribution of ( t , . . . , t n ) (cid:62) as P ( t ≤ a , . . . , t n ≤ a n ) = C ( F x ( a ) , . . . , F x n ( a n )) , where F x i ( t ) = 1 − S x i ( t ) is the cumulative distribution function and the function C is an n -copula used to capture spatial dependence.The current package assumes a spatial version of the Gaussian copula (Li 2010), defined as C ( u , . . . , u n ) = Φ n (cid:0) Φ − { u } , . . . , Φ − { u n } ; R (cid:1) , (12)where Φ n ( · , . . . , · ; R ) denotes the distribution function of N n ( , R ). To allow for a nugget ef-fect, we consider R [ i, j ] = θ ρ ( s i , s j ; θ )+(1 − θ ) I ( s i = s j ), where ρ ( s i , s j ; θ ) = exp {− θ (cid:107) s i − s j (cid:107)} . Here θ ∈ [0 , θ controls the spatial decay overdistance. Note that all the diagonal elements of R are ones, so it is also a correlation matrix.Under the above spatial Gaussian copula, the likelihood function based on upon the completedata { ( t i , x i , s i ) , i = 1 , . . . , n } is L = | R | − / exp (cid:26) − z (cid:62) ( R − − I n ) z (cid:27) n (cid:89) i =1 f x i ( t i ) , where z i = Φ − { F x i ( t i ) } and f x i ( t ) is the density function corresponding to S x i ( t ). We nextdiscuss two marginal spatial survival models for S x i ( t ) that are accommodated in the package.Note that for large n , the FSA introduced in Section 2.1 (with (cid:15) replaced by 1 − θ ) can beapplied. Assume that t i | x i marginally follows the proportional hazards (PH) model with cdf F x i ( t ) = 1 − exp (cid:110) − Λ ( t ) e x (cid:62) i β (cid:111) (13)and density f x i ( t ) = exp (cid:110) − Λ ( t ) e x (cid:62) i β (cid:111) λ ( t ) e x (cid:62) i β , spBayesSurv, version 1.1.3 where β is a p × λ ( t ) is the baseline hazard functionand Λ ( t ) = (cid:82) t λ ( s ) ds is the cumulative baseline hazard function. The piecewise exponentialmodel provides a flexible framework to deal with the baseline hazard (e.g., Walker and Mallick1997). We partition the time period R + into M intervals, say I k = ( d k − , d k ] , k = 1 , . . . , M ,where d = 0 and d M = ∞ . Specifically, we set d k to be the kM th quantile of the empiricaldistribution of the observed survival times for k = 1 , . . . , M −
1. The baseline hazard is thenassumed to be constant within each interval, i.e., λ ( t ) = M (cid:88) k =1 h k I { t ∈ I k } , where h k s are unknown hazard values. Consequently, the cumulative baseline hazard functioncan be written as Λ ( t ) = M ( t ) (cid:88) k =1 h k ∆ k ( t ) , where M ( t ) = min { k : d k ≥ t } and ∆ k ( t ) = min { d k , t } − d k − . After incorporating spatialdependence via the copula in Equation 12, the spCopulaCoxph function considers the followingprior distributions: β ∼ N p ( β , S ) ,h k | h iid ∼ Γ( r h, r ) , k = 1 , . . . , M, ( θ , θ ) ∼ Beta( θ a , θ b ) × Γ( θ a , θ b )The spCopulaCoxph function sets the following default hyperparameter values: M = 10, r = 1, h = ˆ h , β = , S = 10 I p , θ = ( θ a , θ b , θ a , θ b ) (cid:48) = (1 , , , h is themaximum likelihood estimate of the rate parameter from fitting an exponential PH model. Afunction indeptCoxph is also provided to fit the non-spatial standard PH model with abovebaseline and prior settings. The function argument prior allows users to specify these priorparameters in a list with elements defined as follows:element M r0 h0 beta0 S0 theta0 symbol
M r h β S θ We assume that y i = log t i given x i marginally follows a LDDPM model (De Iorio et al. F x i ( t ) = (cid:90) Φ (cid:18) log t − x (cid:62) i β σ (cid:19) dG { β , σ } , (14)where Φ( · ) is the cdf of the standard normal, and G follows the Dirichlet Process (DP) prior.This Bayesian nonparametric model treats the conditional distribution F x as a function-valuedparameter and allows its variance, skewness, modality and other features to flexibly vary withthe x covariates. After incorporating spatial dependence via the copula in Equation 12, the ournal of Statistical Software spCopulaDDP assumes the following prior distributions: G = N (cid:88) k =1 w k δ ( β k ,σ k ) , w k = V k k − (cid:89) j =0 (1 − V j ) , V = 0 , V N = 1 V k iid ∼ Beta(1 , α ) , k = 1 , . . . , N, α ∼ Γ( a , b ) β k | µ iid ∼ N p ( µ , Σ ) , k = 1 , . . . , N, µ ∼ N p ( m , S ) σ − k | Σ iid ∼ Γ( ν a , ν b ) , k = 1 , . . . , N, Σ − ∼ W p (cid:0) ( κ Σ ) − , κ (cid:1) ( θ , θ ) ∼ Beta( θ a , θ b ) × Γ( θ a , θ b ) . The following default hyperparameters are considered in spCopulaDDP : a = b = 2, ν a = 3, ν b = ˆ σ , θ = ( θ a , θ b , θ a , θ b ) (cid:48) = (1 , , , m = ˆ β , S = ˆ Σ , Σ = 30 ˆ Σ , and κ = 7, whereˆ β and ˆ σ are the maximum likelihood estimates of β and σ from fitting the log-normalaccelerated failure time model log( t i ) = x (cid:62) i β + σ(cid:15) i , (cid:15) i ∼ N (0 , Σ is the asymptoticcovariance estimate for ˆ β . A function anovaDDP is also provided to fit the non-spatial LDDPMmodel in Equation 14 with above prior settings. The function argument prior allows usersto specify these prior parameters in a list with elements defined as follows:element N a0 b0 m0 S0 k0 Sig0 theta0 symbol
N a b m S κ Σ θ PH model with spatial copula
The following code is used to fit the piecewise exponential PH model in Equation 13 withthe Gaussian spatial copula in Equation 12 using M = 20 and default priors. We consider K = 100 and B = 1043 for the number of knots and blocks in the FSA of R . The totalrunning time is 15445 seconds. R> set.seed(1)R> mcmc <- list(nburn = 5000, nsave = 2000, nskip = 4, ndisplay = 1000);R> prior <- list(M = 20, nknots = 100, nblock = 1043);R> ptm <- proc.time()R> res1 <- spCopulaCoxph(formula = Surv(time, cens) ~ age + sex + wbc + tpi,+ data = d, mcmc = mcmc, prior = prior,+ Coordinates = cbind(d$xcoord, d$ycoord));R> proc.time() - ptm user system elapsed15262.274 177.716 15444.913
R> (sfit1 <- summary(res1))
Spatial Copula Cox PH model with piecewise constant baseline hazardsCall: spBayesSurv, version 1.1.3 spCopulaCoxph(formula = Surv(time, cens) ~ age + sex + wbc +tpi, data = d, mcmc = mcmc, prior = prior, Coordinates = cbind(d$xcoord,d$ycoord))Posterior inference of regression coefficients(Adaptive M-H acceptance rate: 0.2501):Mean Median Std. Dev. 95%CI-Low 95%CI-Uppage 0.0277864 0.0278065 0.0019297 0.0240332 0.0315580sex 0.0522938 0.0527421 0.0588919 -0.0625843 0.1662136wbc 0.0027808 0.0027899 0.0003767 0.0020071 0.0034546tpi 0.0257918 0.0257969 0.0081385 0.0087972 0.0411955Posterior inference of spatial sill and range parameters(Adaptive M-H acceptance rate: 0.2112):Mean Median Std. Dev. 95%CI-Low 95%CI-Uppsill 0.23051 0.23352 0.05587 0.10222 0.32903range 0.41801 0.34165 0.34272 0.03715 1.31802Log pseudo marginal likelihood: LPML=-5929.357Number of subjects: n=1043 Note that the higher the value of z i = Φ − { F x i ( t i ) } is, the longer the survival time t i (i.e.,lower mortality rate) would be. The posterior sample of z i s is saved in res1$Zpred . Thetrace plots, survival curves, and the map of the posterior mean of z i values can be obtainedusing the code similarly as in Section 2.4. LDDPM model with spatial copula
The following code is used to fit the LDDPM model in Equation 14 with the Gaussian spatialcopula in Equation 12 using N = 10 and default priors. For the FSA, K = 100 and B = 1043are used. The total running time is 20056 seconds. Note there is no summary output as before,as we are fitting a nonparametric model. The trace plots, survival curves, and map of z i s canbe obtained using the same code used for the PH copula model. R> set.seed(1)R> mcmc <- list(nburn = 5000, nsave = 2000, nskip = 4, ndisplay = 1000)R> prior <- list(N = 10, nknots = 100, nblock = 1043)R> ptm <- proc.time()R> res1 <- spCopulaDDP(formula = Surv(time, cens) ~ age + sex + wbc + tpi,+ data = d, mcmc = mcmc, prior = prior,+ Coordinates = cbind(d$xcoord, d$ycoord))R> proc.time() - ptm user system elapsed19876.947 178.595 20056.744
R> sum(log(res1$cpo)); ournal of Statistical Software [1] -5931.5
5. Conclusions
There is a wealth of R packages for non-spatial survival data, starting with survival , includedwith all base installs of R . The survival package fits (discretely) stratified semiparametric PHmodels to right-censored data with exchangeable gamma frailties, as well as left-truncateddata, time-dependent covariates, etc. Parametric log-logistic, Weibull and log-normal AFTmodels can also be fit by this package. From there, there are many packages for variousmodels and types of censoring; a partial review discussing several available R packages isgiven by Zhou and Hanson (2015); also see Zhou and Hanson (2017). In comparison thereare very few R packages for spatially correlated survival data, with the notable exceptions of R2BayesX and spatsurv , both of which focus on PH exclusively. The spBayesSurv packageallows the routine fitting of several popular semiparametric and nonparametric models tospatial survival data. spBayesSurv can also handle non-spatial survival data using either exchangeable Gaussianor no frailty models. Another unintroduced function is survregbayes2 which implementsthe Polya tree based PH, PO, and AFT models of Hanson (2006) and Zhao, Hanson, andCarlin (2009) for areally-referenced data. As pointed out in these papers, MCMC mixingfor Polya tree models can be highly problematic when the true baseline survival function isvery different from the parametric family that centers the Polya tree; the TBP prior providesmuch improved MCMC mixing with essentially the same quality of fit as Polya trees. Anotherfunction very recently added function is
SuperSurvRegBayes , which provides Bayes factorsfor testing among PO, PH, and AFT, as well as three other survival models Zhang, Hanson,and Zhou (2018).Future additions to spBayesSurv include spatial copula (both georeferenced and areal) ver-sions of the PH, PO, and AFT models using TBP priors, as well as continuously-stratifiedproportional hazards and proportional odds models. An extension of all semiparametric mod-els to additive linear structure, which is already incorporated into
BayesX , is also planned.Finally, computational efficiency can be gained by replacing some of the adaptive MCMCupdates with gradient-based updates for the semiparametric models, e.g. the IWLS updatesimplemented in
BayesX for the PH model (Hennerfeind et al.
Acknowledgments
This research is partially funded by Grant R03CA176739 from National Institutes of Health.The authors would like to thank referees for their valuable comments, and all users who havereported bugs and given suggestions.
References
Antoniak CE (1974). “Mixtures of Dirichlet Processes With Applications to Bayesian Non-parametric Problems.”
The Annals of Statistics , , 1152–1174.8 spBayesSurv, version 1.1.3 Arbia G, Espa G, Giuliani D, Micciolo (2016). “A Spatial Analysis of Health and Pharma-ceutical Firm Survival.”
Journal of Applied Statistics , p. in press.Banerjee S, Carlin BP, Gelfand AE (2014).
Hierarchical Modeling and Analysis for SpatialData, Second Edition . Chapman and Hall/CRC Press.Banerjee S, Dey DK (2005). “Semiparametric Proportional Odds Models for Spatially Corre-lated Survival Data.”
Lifetime Data Analysis , (2), 175–191.Banerjee S, Gelfand AE, Finley AO, Sang H (2008). “Gaussian Predictive Process Models forLarge Spatial Data Sets.” Journal of the Royal Statistical Society B , (4), 825–848.B´ardossy A (2006). “Copula-Based Geostatistical Models for Groundwater Quality Parame-ters.” Water Resources Research , (11), 1–12.Belitz C, Brezger A, Klein N, Kneib T, Lang S, Umlauf N (2015). BayesX - Software forBayesian Inference in Structured Additive Regression Models
Journalof the Royal Statistical Society B , (2), 192–236.Carlin BP, Louis TA (2010). Bayes and Empirical Bayes Methods for Data Analysis . Chapmanand Hall/CRC.Chen Y, Hanson T, Zhang J (2014). “Accelerated Hazards Model Based on ParametricFamilies Generalized With Bernstein Polynomials.”
Biometrics , (1), 192–201.Chiou SH, Kang S, Yan J (2015). “Semiparametric Accelerated Failure Time Modeling forClustered Failure Times From Stratified Sampling.” Journal of the American StatisticalAssociation , (510), 621–629.Cox DR, Oakes D (1984). Analysis of Survival Data . Chapman & Hall: London.Cox DR, Snell EJ (1968). “A General Definition of Residuals.”
Journal of the Royal StatisticalSociety B , (2), 248–275.Darmofal D (2009). “Bayesian Spatial Survival Models for Political Event Processes.” Amer-ican Journal of Political Science , (1), 241–257. ISSN 1540-5907.De Iorio M, Johnson WO, M¨uller P, Rosner GL (2009). “Bayesian Nonparametric Nonpro-portional Hazards Survival Modeling.” Biometrics , (3), 762–771.Finley AO, Sang H, Banerjee S, Gelfand AE (2009). “Improving the Performance of PredictiveProcess Modeling for Large Datasets.” Computational Statistics & Data Analysis , (8),2873–2884.Geisser S, Eddy WF (1979). “A Predictive Approach to Model Selection.” Journal of theAmerican Statistical Association , (365), 153–160.Gelman A (2006). “Prior Distributions for Variance Parameters in Hierarchical Models (Com-ment on Article by Browne and Draper).” Bayesian analysis , (3), 515–534. ournal of Statistical Software Bernoulli , (2), 223–242.Hanson T, Johnson W, Laud P (2009). “Semiparametric Inference for Survival Models WithStep Process Covariates.” Canadian Journal of Statistics , (1), 60–79.Hanson TE (2006). “Inference for Mixtures of Finite Polya Tree Models.” Journal of theAmerican Statistical Association , (476), 1548–1565.Hanson TE, Branscum AJ, Johnson WO (2014). “Informative g -Priors for Logistic Regres-sion.” Bayesian Analysis , (3), 597–612.Henderson R, Shimakura S, Gorst D (2002). “Modeling Spatial Variation in Leukemia SurvivalData.” Journal of the American Statistical Association , (460), 965–972.Hennerfeind A, Brezger A, Fahrmeir L (2006). “Geoadditive Survival Models.” Journal of theAmerican Statistical Association , (475), 1065–1075.Johnson ME, Moore LM, Ylvisaker D (1990). “Minimax and Maximin Distance Designs.” Journal of statistical planning and inference , (2), 131–148.Kammann EE, Wand MP (2003). “Geoadditive Models.” Applied Statistics , , 1–18.Kneib T (2006). “Mixed Model-Based Inference in Geoadditive Hazard Regression for Interval-Censored Survival Times.” Computational Statistics & Data Analysis , (2), 777–792.Kneib T, Fahrmeir L (2007). “A Mixed Model Approach for Geoadditive Hazard Regression.” Scandinavian Journal of Statistics , (1), 207–228.Konomi BA, Sang H, Mallick BK (2014). “Adaptive Bayesian Nonstationary Modeling forLarge Spatial Datasets Using Covariance Approximations.” Journal of Computational andGraphical Statistics , , 802–929.Kuo L, Mallick B (1998). “Variable Selection for Regression Models.” Sankhy¯a: The IndianJournal of Statistics, Series B , , 65–81.Lavine M (1992). “Some Aspects of Polya Tree Distributions for Statistical Modelling.” TheAnnals of Statistics , , 1222–1235.Lavine ML, Hodges JS (2012). “On Rigorous Specification of ICAR Models.” The AmericanStatistician , (1), 42–49.Li J (2010). Application of Copulas as a New Geostatistical Tool . Ph.D. thesis, Institut f¨urWasser- und Umweltsystemmodellierung.Li J, Hong Y, Thapa R, Burkhart HE (2015). “Survival Analysis of Loblolly Pine Trees WithSpatially Correlated Random Effects.”
Journal of the American Statistical Association , (510), 486–502.Li Y, Lin X (2006). “Semiparametric Normal Transformation Models for Spatially CorrelatedSurvival Data.” Journal of the American Statistical Association , (474), 591–603.0 spBayesSurv, version 1.1.3 M¨uller P, Quintana F, Jara A, Hanson T (2015).
Bayesian Nonparametric Data Analysis .Springer-Verlag: New York.Nardi A, Schemper M (2003). “Comparing Cox and Parametric Models in Clinical Studies.”
Statistics in Medicine , (23), 3597–3610.Nychka D, Furrer R, Paige J, Sain S (2015). “fields: Tools for Spatial Data.” doi:10.5065/D6W957CT . R package version 8.10, URL .Sang H, Huang JZ (2012). “A Full Scale Approximation of Covariance Functions for LargeSpatial Data Sets.” Journal of the Royal Statistical Society B , (1), 111–132.Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A (2002). “Bayesian Measures of ModelComplexity and Fit.” Journal of the Royal Statistical Society B , (4), 583–639.Taylor BM (2015). “Auxiliary Variable Markov Chain Monte Carlo for Spatial Survival andGeostatistical Models.” arXiv preprint arXiv:1501.01665 .Taylor BM (2017). “Spatial Modelling of Emergency Service Response Times.” Journal of theRoyal Statistical Society A , (2), 433–453.Taylor BM, Rowlingson BS (2017). “Spatsurv: An R Package for Bayesian Inference WithSpatial Survival Models.” Journal of Statistical Software , (4), 1–32.Therneau T, Crowson C, Atkinson E (2017). Using Time Dependent Covariates and TimeDependent Coefficients in the Cox Model . URL http://cran.es.r-project.org/web/packages/survival/vignettes/timedep.pdf .Therneau TM (2015).
A Package for Survival Analysis in S . Version 2.38, URL https://CRAN.R-project.org/package=survival .Turnbull BW (1974). “Nonparametric Estimation of a Survivorship Function With DoublyCensored Data.”
Journal of the American Statistical Association , (345), 169–173.Umlauf N, Adler D, Kneib T, Lang S, Zeileis A (2015). “Structured Additive RegressionModels: An R Interface to BayesX.” Journal of Statistical Software , (21), 1–46.Verdinelli I, Wasserman L (1995). “Computing Bayes Factors Using a Generalization of theSavage-Dickey Density Ratio.” Journal of the American Statistical Association , (430),614–618.Walker SG, Mallick BK (1997). “Hierarchical Generalized Linear Models and Frailty ModelsWith Bayesian Nonparametric Mixing.” Journal of the Royal Statistical Society B , ,845–860.Waller LA, Gotway CA (2004). Applied Spatial Statistics for Public Health Data . John Wiley& Sons.Wang S, Zhang J, Lawson AB (2012). “A Bayesian Normal Mixture Accelerated Failure TimeSpatial Model and Its Application to Prostate Cancer.”
Statistical Methods in MedicalResearch , http://dx.doi.org/10.1177/0962280212466189 . ournal of Statistical Software Journal of Machine LearningResearch , (Dec), 3571–3594.Wood S (2017). Generalized Additive Models: An Introduction with R . 2 edition. Chapmanand Hall/CRC.Zhang J, Hanson T, Zhou H (2018). “Bayes Factors for Choosing Among Six Common SurvivalModels.”
Lifetime Data Analysis , pp. 1–19. doi:10.1007/s10985-018-9429-4 .Zhao L, Hanson TE, Carlin BP (2009). “Mixtures of Polya Trees for Flexible Spatial FrailtySurvival Modelling.”
Biometrika , (2), 263–276.Zhou H, Hanson T (2015). “Bayesian Spatial Survival Models.” In Nonparametric BayesianInference in Biostatistics , pp. 215–246. Springer-Verlag.Zhou H, Hanson T (2017). “A Unified Framework for Fitting Bayesian Semiparametric Modelsto Arbitrarily Censored Survival Data, Including Spatially-Referenced Data.”
Journal ofthe American Statistical Association , in press .Zhou H, Hanson T (2018). spBayesSurv: Bayesian Modeling and Analysis of Spatially Cor-related Survival Data . R package version > = 1.1.3, URL https://CRAN.R-project.org/package=spBayesSurv .Zhou H, Hanson T, Jara A, Zhang J (2015a). “Modeling County Level Breast Cancer SurvivalData Using a Covariate-Adjusted Frailty Proportional Hazards Model.” The Annals ofApplied Statistics , (1), 43–68.Zhou H, Hanson T, Knapp R (2015b). “Marginal Bayesian Nonparametric Model for Timeto Disease Arrival of Threatened Amphibian Populations.” Biometrics , (4), 1101–1110.Zhou H, Hanson T, Zhang J (2017). “Generalized Accelerated Failure Time Spatial FrailtyModel for Arbitrarily Censored Data.” Lifetime Data Analysis , (3), 495–515. Affiliation:
Haiming ZhouDivision of StatisticsNorthern Illinois UniversityE-mail: [email protected]
Timothy HansonStrategic & Scientific OperationsMedtronic Inc., Minneapolis, Minnesota, U.S.A.E-mail: [email protected] spBayesSurv, version 1.1.3 Jiajia ZhangDepartment of Epidemiology and BiostatisticsUniversity of South CarolinaE-mail: [email protected]
Journal of Statistical Software published by the Foundation for Open Access Statistics
MMMMMM YYYY, Volume VV, Issue II
Submitted: yyyy-mm-dd doi:10.18637/jss.v000.i00