DPpackage (version 1.1-7.4)

LDPDdoublyint: Bayesian analysis for a Linear Dependent Poisson Dirichlet Process Mixture Models for the Analysis of Doubly-Interval-Censored Data

Description

This function generates a posterior density sample for a linear dependent Poisson Dirichlet process mixture model for the analysis of doubly-invertval-censored survival data.

Usage

LDPDdoublyint(onset,failure,p,xonset,q,xfailure,xpred,
              grid,prior,mcmc,state,status,work.dir=NULL)

Arguments

onset

a matrix giving the interval limits for the onset times. For multiple variables the limits must be included in a sequential manner for each variable, i.e., (L1,U1,L2,U2,...).

failure

a matrix giving the interval limits for the time-to-event times. For multiple variables the limits must be included in a sequential manner for each variable, i.e. (L1,U1,L2,U2,...).

p

an integer giving the number of predictors included for each onset variable. Different design matrices are allowed for the different responses but of the same p-dimension.

xonset

a matrix giving the design matrices for each onset time. For multiple variables the design matrices must be included in order and includes the intercepts, i.e. (XO1,XO2,...).

q

an integer giving the number of predictors included for each time-to-event variable. Different design matrices are allowed for the different responses but of the same q-dimension.

xfailure

a matrix giving the design matrices for each time-to-event variable. For multiple variables the design matrices must be included in order and includes the intercepts, i.e. (XT1,XT2,...).

xpred

a matrix giving the value of the predictors for which survival and hazard functions will be evaluated. The number of columns of xpred should be (p+q)*nvar/2 where nvar is the number of variables. The design matrices for the predictions must include onset predictors first and then time-to-event predictors, i.e., (XO1,XO2,...,XT1,XT2,...).

grid

a matrix including the grids where survival and hazard functions are evaluated. Each row must include the grid points for different variable.

prior

a list giving the prior information. The list includes the following parameter: q, a0 and b0 giving the hyperparameters for prior distribution of the a precision parameter of the Poisson-Dirichlet process prior, mub and sigmab giving the hyperparameters for the prior distributions of the b precision parameter of the Poisson-Dirichlet process prior, nu and tinv giving the hyperparameters of the inverted Wishart prior distribution for the kernel covariance matrix, mb and Sb giving the hyperparameters of the normal prior distribution for the mean of the normal baseline distribution, nub and tbinv giving the hyperparameters of the inverted Wishart prior distribution for the for the covariance matrix of the normal baseline distribution, and maxm giving the finite truncation limit of the sitck-breaking representation of the Poisson-Dirichlet process.

mcmc

a list giving the MCMC parameters. The list must include the following integers: nburn giving the number of burn-in scans, nskip giving the thinning interval, nsave giving the total number of scans to be saved, ndisplay giving the number of saved scans to be displayed on screen (the function reports on the screen when every ndisplay iterations have been carried out), and tune1 and tune2 the parameters needed for the MH candidate generating distribution for the precision parameters of the Poisson-Dirichlet process prior.

state

a list giving the current value of the parameters. This list is used if the current analysis is the continuation of a previous analysis.

status

a logical variable indicating whether this run is new (TRUE) or the continuation of a previous analysis (FALSE). In the latter case the current value of the parameters must be specified in the object state.

work.dir

working directory.

Value

An object of class LDPDdoublyint representing the LDPD mixture of survival models fit. Generic functions such as print, plot, summary, and predict have methods to show the results of the fit. The results include mb, Sb, sigma, the precision parameter alpha, and the number of clusters.

The list state in the output object contains the current value of the parameters necessary to restart the analysis. If you want to specify different starting values to run multiple chains set status=TRUE and create the list state based on this starting values. In this case the list state must include the following objects:

alpha

a vector giving the value of the precision parameters a and b.

b

a matrix of dimension maxm times the number of columns in the design matrix ((p+q)*nvar/2), giving the regression coefficients for each cluster.

sigma

a matrix of dimension nvar times nvar giving the kernel covariance matrix.

mb

giving the mean of the normal baseline distributions.

Sb

giving the covariance matrix the normal baseline distributions.

ncluster

an integer giving the number of clusters.

ss

an interger vector defining to which of the ncluster clusters each subject belongs.

y

a matrix of dimension nsubject times nvar giving the value of the imputed log-survival times.

Details

This generic function fits a Linear Dependent Poisson-Dirichlet Process Mixture of Survival models as in Jara et al. (2010). The joint distribution of the true chronological onset times and true time-to-events, \(T_i\), is modeled using the mixture model $$log(T_i)=y_i | f_{X_i} \sim f_{X_i}$$ $$f_{X_i} = \int N(X_i \beta, \Sigma) dG(\beta)$$ $$G | a, b, G_0 \sim PD(a,b,G_0)$$ where, \(G_0 = N(\beta| m_b, S_b)\). To complete the model specification, independent hyperpriors are assumed, $$\Sigma | \nu, T \sim IW(\nu,T)$$ $$a | q, a_0, b_0 \sim q \delta_0 + (1-q)Beta(a_0,b_0)$$ $$b | a, \mu_b, \sigma_b \sim N(\mu_b,\sigma_b) I(-a,\infty)$$ $$m_b | m_0, S_0 \sim N(m_0,S_0)$$ $$S_b | \nu_b, Tb \sim IW(\nu_b,T_b)$$

Note that the inverted-Wishart prior is parametrized such that if \(A \sim IW_q(\nu, \psi)\) then \(E(A)= \psi^{-1}/(\nu-q-1)\).

The precision parameters are updated using a MH step. The candidate generating distribution is of the form $$a' | a, tune_1 \sim 0.5 \delta_0 + 0.5 N(a,tune_1^2)$$ $$b' | b, a', tune_2 \sim N(b,tune_2^2) I(-a',\infty)$$

The computational implementation of the model is based on the maxm truncation of stick-breaking representation of the model (see, Jara et al., 2009).

References

Jara, A., Lesaffre, E., De Iorio, M., Quintana, F. (2010). Bayesian semiparametric inference for multivariate doubly-interval-censored data. Annals of Applied Statistics, 4: 2126-2149.

See Also

LDDPdensity, LDDPsurvival

Examples

# NOT RUN {
###############################################################
# HIV-AIDS Data
###############################################################

  data(hiv)
  attach(hiv)

###############################################################
# Working folder
###############################################################

  work.dir <- NULL

###############################################################
# Response and design matrices
###############################################################

  nsubject <- length(onsetL)
  onset <- cbind(onsetL,onsetU)
  failure <- cbind(failureL,failureU)

  intercept <- rep(1,nsubject)
  p <- 2
  xonset <- cbind(IntO=intercept,trtO=trt)
  q <- 2
  xfailure <- cbind(IntF=intercept,trtF=trt)

###############################################################
# Predictions
###############################################################

   grid <- matrix(c(rep(seq(0,30,len=30),1),
                    rep(seq(0,20,len=30),1)),nrow=2,byrow=T)

   xpred <- matrix(c(rep(c(1,0),1),rep(c(1,0),1),
                     rep(c(1,1),1),rep(c(1,1),1)),
                     nrow=2,byrow=T)

   colnames(xpred) <- colnames(cbind(xonset,xfailure))

###############################################################
# Initial state
###############################################################

  state <- NULL

###############################################################
# Prior
###############################################################
  
  prior<-list(a0=1,
              b0=1,
              q=0.5,
              mub=10,
              sigmab=200,
              nu=4,
              tinv=diag(1,2),
              nub=6,
              tbinv=diag(1,4),
              m0=rep(0,4),
              S0=diag(100,4),
              maxm=40)

###############################################################
# MCMC
###############################################################
  
  mcmc<-list(nburn=5000,
             nskip=9,
             ndisplay=100,
             nsave=5000,
             tune1=0.25,
             tune2=1)

###############################################################
# Fitting the Model
###############################################################

  fit1 <- LDPDdoublyint(onset=onset,failure=failure,
                        p=p,xonset=xonset,
                        q=q,xfailure=xfailure,
                        xpred=xpred,grid=grid,
                        prior=prior,
                        mcmc=mcmc,
                        state=state,
                        status=TRUE,
                        work.dir=work.dir)
 
  fit1
  summary(fit1)

###############################################################
# Getting Information for Predictions
###############################################################
  
# Without CI bands and intervals

  fit1.pred <- predict(fit1)
  fit1.pred
  plot(fit1.pred)

# With CI bands and intervals
 
  fit1.pred <- predict(fit1,compute.band=TRUE)
  fit1.pred
  plot(fit1.pred)
# }