LDPDdoublyint: Bayesian analysis for a Linear Dependent Poisson Dirichlet...
In DPpackage: Bayesian Nonparametric Modeling in R

Description Usage Arguments Details Value Author(s) References See Also Examples

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.

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

`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.

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, Ti, is modeled using the mixture model

log(Ti)=yi | fXi ~ fXi

fXi = \int N(Xi beta, Sigma) dG(beta)

G | a, b, G0 ~ PD(a, b, G0)

where, G0 = N(beta | mb, Sb). To complete the model specification, independent hyperpriors are assumed,

Sigma | nu, T ~ IW(nu,T)

a | q, a0, b0 ~ q delta0 + (1-q) Beta(a0,b0)

b | a, mub, sigmab ~ N(mub,sigmab)I(-a,infty)

mb | m0, S0 ~ N(m0,S0)

Sb | nub, Tb ~ IW(nub,Tb)

Note that the inverted-Wishart prior is parametrized such that if A ~ IWq(nu, psi) then E(A)= psiinv/(nu-q-1).

The precision parameters are updated using a MH step. The candidate generating distribution is of the form

a' | a, tune1 ~ 0.5 delta0 + 0.5 N(a,tune1**2)

b' | b, a', tune2 ~ N(b,tune2**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).

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.

Alejandro Jara <atjara@uc.cl>

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.

LDDPdensity, LDDPsurvival

## 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)

## End(Not run)